Industrial robot motion trail obstacle avoidance method based on deep learning

Abstract

The invention discloses an industrial robot motion trail obstacle avoidance method based on deep learning, which relates to the technical field of industrial robots and comprises the following steps: extracting a position sequence, a speed sequence and an acceleration sequence from trail data of an industrial robot motion trail; respectively generating a geometric shape feature vector and a dynamic feature vector through a convolutional neural network, and fusing the geometric shape feature vector and the dynamic feature vector to obtain an initial representation vector; according to the industrial robot motion trail obstacle avoidance method based on deep learning, it is ensured that the finally generated industrial robot motion trail is kept coordinated and unified among the geometric structure, the dynamic smoothness and the obstacle avoidance reliability, and therefore the operation stability, the operation efficiency and the safety level of an industrial robot in a complex production environment are improved.

Description

Technical Field

[0001] This invention relates to the field of industrial robot technology, and specifically to a method for obstacle avoidance of industrial robot motion trajectory based on deep learning. Background Technology

[0002] With the rapid development of intelligent manufacturing and flexible production lines, industrial robots are increasingly being used in complex work environments such as welding, assembly, handling, and grinding. In industrial settings characterized by multi-device collaboration, densely packed workstations, and frequent dynamic obstacles, robots need to plan safe, smooth, and efficient motion trajectories in real time to avoid collisions with workpieces, fixtures, or other equipment. Therefore, achieving reliable obstacle avoidance while meeting kinematic and dynamic constraints has become a key technical challenge for improving the autonomous operation capabilities and production efficiency of industrial robots.

[0003] Most existing obstacle avoidance methods for industrial robots are based on geometric space modeling and collision detection mechanisms. They achieve path correction by calculating the minimum distance between the robot and obstacles or constructing a safe envelope region; or they introduce simple kinematic constraints such as velocity and acceleration to improve trajectory smoothness based on traditional trajectory planning algorithms. However, these methods typically focus on obstacle avoidance judgment at the geometric space level, or only optimize dynamic indicators individually, lacking a systematic analysis of the coupling relationship between trajectory geometry and dynamic behavior. This makes it difficult to comprehensively characterize the real motion characteristics of robots in complex industrial environments. From a data structure perspective, the motion trajectory of an industrial robot is essentially high-dimensional time-series data, containing multi-order derivative information such as position, velocity, and acceleration. The curvature changes and torsion distribution of the trajectory in space reflect its geometric shape characteristics, while the acceleration change trend and its temporal evolution reflect the dynamic execution state. In actual motion, there is often a significant coupling relationship between geometric and dynamic characteristics: abrupt curvature changes may cause acceleration peaks or abnormal acceleration rates of change, while dynamic constraints, in turn, limit the range of achievable geometric paths. If these two types of information cannot be effectively integrated during trajectory representation, it can easily lead to trajectory distortion, resulting in theoretically feasible planning results but causing problems such as vibration, impact, or excessive energy consumption during actual execution. Furthermore, existing methods, when processing high-dimensional trajectory data, often employ fixed feature extraction or manual empirical parameter modeling, lacking the ability to automatically learn the deep structural features of the trajectory, making it difficult to achieve adaptive optimization in complex environments. Simultaneously, during trajectory adjustment and obstacle avoidance decisions, there is often a lack of feedback mechanisms for the overall structural consistency of the trajectory. After multiple local corrections, global structural imbalances or disruptions in dynamic continuity may occur, affecting the robot's operational stability and safety. Summary of the Invention

[0004] The purpose of this invention is to provide a method for obstacle avoidance of industrial robot motion trajectory based on deep learning, thereby solving the problems existing in the prior art.

[0005] To achieve the above objectives, the present invention provides the following technical solution: a method for obstacle avoidance of industrial robot motion trajectory based on deep learning, comprising: S1, extracting position sequence, velocity sequence, and acceleration sequence from the trajectory data of the industrial robot motion trajectory, generating geometric feature vector and dynamic feature vector respectively through a convolutional neural network, and fusing the two to obtain an initial representation vector; S2, calculating the synchronization degree between local curvature change and acceleration peak, the spatial overlap density between curvature extreme points and velocity change points, and the directional consistency between trajectory torsion gradient and acceleration rate of change vector for the initial representation vector; if the weighted average of the three indicators exceeds a preset threshold, then the geometric feature vector is determined to be obstructed. There is a strong coupling relationship between shape and dynamic characteristics; S3, based on the initial representation vector determined to be strongly coupled, principal component analysis is used to extract the geometric dynamic coupling principal direction that retains the maximum variance, the joint principal component axis in the cross-modal information set, and the subspace that minimizes the joint reconstruction error of geometric distortion and dynamic distortion after projection, to obtain the dimension-reduced joint feature representation; S4, the sequence of industrial robot trajectory sample points is restored from the dimension-reduced joint feature representation, and input into the long short-term memory network to capture the temporal pattern corresponding to the matching degree of curvature zero-crossing frequency and velocity zero-crossing frequency, as well as the matching trend of the distribution of bending segment length and the statistical distribution of acceleration amplitude, to obtain the trajectory dynamic mode vector.

[0006] Preferably, step S1 includes acquiring trajectory data during the operation of the industrial robot, parsing the trajectory data to obtain a position sequence, a velocity sequence, and an acceleration sequence; inputting the position sequence into a spatial feature extraction convolutional neural network to generate a geometric shape feature vector, and inputting the velocity sequence and acceleration sequence into a dynamic feature extraction convolutional neural network to generate a dynamic feature vector; performing dimensional mapping and feature concatenation on the geometric shape feature vector and the dynamic feature vector, and fusing them to obtain an initial representation vector containing complete spatiotemporal information.

[0007] Preferably, step S2 includes parsing the initial representation vector to obtain geometric and dynamic feature components; extracting local curvature change sequences and acceleration peak sequences from the geometric and dynamic feature components; calculating the synchronization degree values ​​of the local curvature change sequences and acceleration peak sequences; determining the spatial overlap density of curvature extreme points and velocity abrupt change points based on the synchronization degree values; calculating the directional consistency index based on the spatial overlap density; and determining that a strong coupling relationship exists between the geometric shape and dynamic characteristics if the comprehensive coupling coefficient generated by the synchronization degree values, spatial overlap density, and directional consistency index meets the conditions.

[0008] Preferably, step S3 includes obtaining the initial representation vector selected by strong coupling determination, constructing a joint covariance matrix using the initial representation vector; decomposing the joint covariance matrix to extract the geometric dynamic coupling principal direction that retains the maximum variance, and identifying the joint principal component axis in the cross-modal information set; constructing an initial subspace based on the geometric dynamic coupling principal direction and the joint principal component axis, and optimizing the target subspace by minimizing the joint reconstruction error function; mapping the initial representation vector to the target subspace to obtain the dimensionality-reduced joint feature representation.

[0009] Preferably, step S4 includes obtaining the dimensionality-reduced joint feature representation, reconstructing the industrial robot trajectory sample point sequence through inverse projection; calculating the frequency matching degree of curvature and velocity zero-crossing frequency, and the distribution matching trend of bending segment length and acceleration amplitude based on the industrial robot trajectory sample point sequence; combining the frequency matching degree and distribution matching trend to generate a temporal feature sequence, and inputting the temporal feature sequence into a long short-term memory network; obtaining the hidden layer state vector output by the long short-term memory network, and mapping the hidden layer state vector to obtain the trajectory dynamics mode vector.

[0010] Preferably, the method further includes S5: calculating the correspondence score between the path fractal dimension and the dynamic energy dissipation rate for the trajectory dynamics mode vector and the geometric shape feature vector; if the correspondence score is lower than a preset threshold, adjusting the positions of key control points in the industrial robot trajectory sample point sequence through gradient descent iteratively until the adjusted sequence satisfies the predefined kinematic continuity constraints; specifically, obtaining the trajectory dynamics mode vector and the geometric shape feature vector, calculating the path fractal dimension and the dynamic energy dissipation rate respectively; calculating the mutual information value of the path fractal dimension and the dynamic energy dissipation rate as the correspondence score; if the correspondence score is lower than the preset threshold, updating the control point positions along the gradient vector that maximizes the correspondence score to reconstruct the corrected industrial robot trajectory sample point sequence; performing a third-order differential on the corrected industrial robot trajectory sample point sequence to obtain the jerk sequence; if the jerk sequence satisfies the smoothness limit value, determining that the adjusted sequence satisfies the predefined kinematic continuity constraints.

[0011] Preferably, the method further includes S6: based on the adjusted industrial robot trajectory sequence that satisfies kinematic continuity constraints, and combined with the obstacle model within the industrial robot's workspace, calculating the interference probability caused by the synchronization degree between local curvature changes and acceleration peaks within the potential obstacle region; if the interference probability is higher than a preset threshold, then using a dynamic programming algorithm to search for candidate node sequences that avoid the interference probability region, thereby obtaining the optimal obstacle avoidance path. Specifically, this includes obtaining the industrial robot trajectory sequence and obstacle model that satisfy kinematic continuity constraints, identifying potential obstacle regions located near the obstacle model in the trajectory sequence; calculating the synchronization degree based on local curvature changes and acceleration peaks within the potential obstacle region to determine the interference probability; if the interference probability is higher than a preset threshold, then using a dynamic programming algorithm to search for the optimal state transition path passing through candidate nodes outside the potential obstacle region, thereby obtaining the optimal obstacle avoidance path.

[0012] Preferably, the process also includes S7: extracting key node positions with high spatial overlap density between curvature extrema and velocity mutation points from the obstacle avoidance optimization path; verifying the integrity of geometric self-similarity and acceleration sequence autocorrelation peak values ​​using a convolutional neural network to determine the final industrial robot motion trajectory representation; specifically, calculating the spatial overlap density between curvature extrema and velocity mutation points in the obstacle avoidance optimization path; extracting local trajectory segments based on the spatial overlap density; obtaining the geometric self-similarity and acceleration sequence autocorrelation peak values ​​of the local trajectory segments; constructing a multidimensional feature tensor based on the geometric self-similarity and acceleration sequence autocorrelation peak values; performing integrity verification on the multidimensional feature tensor using a convolutional neural network and outputting the integrity verification result; and performing smooth reconstruction on the local trajectory segments based on the integrity verification result to determine the final industrial robot motion trajectory representation.

[0013] Preferably, the process also includes S8, calculating the consistency score between the final industrial robot motion trajectory representation and the original high-dimensional trajectory data of the industrial robot, specifically the alignment degree between the geometric symmetry axis and the velocity vector flipping moment, and the matching degree between the curvature zero-crossing frequency and the velocity zero-crossing frequency. If the consistency score is lower than a preset threshold, the joint principal component axis weights in the initial representation vector are updated through error feedback, and the process re-enters the dimensionality reduction process. Specifically, this includes extracting the geometric symmetry axis of the final industrial robot motion trajectory representation and the velocity vector flipping moment of the original trajectory data to calculate the axis alignment degree; and extracting the curvature zero-crossing frequency of the trajectory representation and the velocity zero-crossing frequency of the original trajectory data to calculate the frequency matching degree.

[0014] Preferably, step S8 further includes calculating a consistency score based on the axis alignment degree and frequency matching degree; if the consistency score is lower than a preset threshold, an error feedback signal is generated to update the joint principal component axis weights in the initial representation vector, and the vector re-enters the dimensionality reduction process based on the updated joint principal component axis weights.

[0015] As can be seen from the above technical solution, the present invention has the following beneficial effects:

[0016] This deep learning-based obstacle avoidance method for industrial robots utilizes convolutional neural networks for adaptive extraction of multidimensional features and a coupled decision mechanism to identify the strong correlation between geometric and dynamic characteristics. Dimensionality reduction is achieved within a unified representation framework, effectively avoiding trajectory distortion caused by the fragmented modeling of geometric and dynamic information in traditional methods. Furthermore, by incorporating a long short-term memory network to learn the temporal patterns of the trajectory, it accurately characterizes the dynamic evolution relationship between curvature and acceleration changes, improving the ability to recognize complex motion patterns. Finally, by introducing an evaluation mechanism that establishes a correspondence between path fractal dimension and dynamic energy dissipation rate, the rationality of the trajectory structure is determined. The system determines the parameters and optimizes key control points using gradient descent when preset conditions are not met, ensuring the trajectory satisfies kinematic continuity constraints. During obstacle avoidance, it constructs an obstacle model of the industrial robot's workspace and calculates the interference probability, then uses dynamic programming to search for low-interference-risk paths, improving the accuracy and safety of obstacle avoidance decisions. Furthermore, it iteratively updates the joint principal component axis weights through a consistency score feedback mechanism, forming a closed-loop optimization structure. This ensures that the final generated industrial robot trajectory maintains coordination and unity between geometric structure, dynamic smoothness, and obstacle avoidance reliability, thereby improving the operational stability, efficiency, and safety of the industrial robot in complex production environments. Attached Figure Description

[0017] Figure 1 This is a flowchart of the obstacle avoidance method for the motion trajectory of an industrial robot according to the present invention. Detailed Implementation

[0018] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0019] like Figure 1As shown, this invention provides a technical solution: a deep learning-based obstacle avoidance method for industrial robot motion trajectories, comprising: S1, extracting position sequence, velocity sequence, and acceleration sequence from the trajectory data of the industrial robot motion trajectory, generating geometric feature vectors and dynamic feature vectors respectively through a convolutional neural network, and fusing the two to obtain an initial representation vector; S2, calculating the synchronization degree between local curvature changes and acceleration peaks, the spatial overlap density between curvature extreme points and velocity abrupt change points, and the directional consistency between the trajectory torsion gradient and the acceleration rate of change vector for the initial representation vector; if the weighted average of the three indicators exceeds a preset threshold, it is determined that there is a strong coupling between the geometric shape and dynamic characteristics. S3. Based on the initial representation vector determined to be strongly coupled, principal component analysis is used to extract the geometric-dynamic coupling principal direction that retains the maximum variance, the joint principal component axis in the cross-modal information set, and the low-dimensional subspace that minimizes the joint reconstruction error of geometric and dynamic distortion after projection, to obtain the dimensionality-reduced joint feature representation; S4. The trajectory sample point sequence of the industrial robot is reconstructed from the dimensionality-reduced joint feature representation, and input into a long short-term memory network to capture the temporal pattern corresponding to the matching degree of curvature zero-crossing frequency and velocity zero-crossing frequency, as well as the matching trend of the distribution of bending segment length and the statistical distribution of acceleration amplitude, to obtain the trajectory dynamics mode vector; S5. For the trajectory dynamics mode vector and the geometric-dynamic coupling principal direction that retains the maximum variance, principal component analysis is used to extract the joint principal component axis in the cross-modal information set, and the low-dimensional subspace that minimizes the joint reconstruction error of geometric and dynamic distortion after projection, to obtain the trajectory dynamics mode vector; The path fractal dimension and dynamic energy dissipation rate are correlated using the shape feature vector. If the correlation score is lower than a preset threshold, the positions of key control points in the industrial robot trajectory sample point sequence are iteratively adjusted using gradient descent until the adjusted sequence satisfies the predefined kinematic continuity constraints. S6. Based on the adjusted industrial robot trajectory sequence that satisfies the kinematic continuity constraints, and combined with the obstacle model within the industrial robot's workspace, the interference probability caused by the synchronization degree of local curvature changes and acceleration peaks within the potential obstacle region is calculated. If the interference probability is higher than a preset threshold, a dynamic programming algorithm is used to search for candidate node sequences that avoid the interference probability region. S7. Extract key node positions with high spatial overlap density between curvature extrema and velocity mutation points from the obstacle avoidance optimization path. Verify the integrity of geometric self-similarity and acceleration sequence autocorrelation peaks through convolutional neural networks to determine the final industrial robot motion trajectory representation. S8. Calculate the consistency score between the final industrial robot motion trajectory representation and the original high-dimensional trajectory data of the industrial robot, based on the alignment degree of the geometric symmetry axis and the velocity vector flipping moment, as well as the matching degree of curvature zero-crossing frequency and velocity zero-crossing frequency. If the consistency score is lower than the preset threshold, update the joint principal component axis weights in the initial representation vector through error feedback and re-enter the dimensionality reduction process.

[0020] In the above scheme, the raw trajectory data of the industrial robot's end effector is first uniformly sampled and time-aligned. Continuous position points are acquired at fixed time intervals, and the corresponding velocity and acceleration change sequences are calculated based on the changing trends of adjacent sampling points. In this way, the original discrete trajectory is transformed into multi-dimensional time-series data containing spatial shape information and motion state information.

[0021] Subsequently, two feature extraction channels are constructed. The first channel takes the position sequence as input and extracts local curvature changes, spatial turning distribution, and overall path morphology features through a multi-layer convolutional structure, forming a geometric shape feature vector. The second channel takes the velocity and acceleration sequences as input and extracts motion abrupt change features, acceleration peak distribution, and dynamic change trend features, forming a dynamic feature vector. Through feature concatenation and linear mapping, the two types of features are fused into an initial representation vector of a unified dimension, thereby achieving a preliminary fusion expression of geometric and dynamic information.

[0022] Based on the fusion expression, a multi-dimensional coupling evaluation mechanism is constructed. By analyzing the synchronization of local bending changes and acceleration peaks in the time dimension, the degree of overlap between bending extreme points and velocity change points in space is statistically analyzed, and the consistency between the trajectory spatial change trend and the acceleration change direction is compared, comprehensively evaluating the coupling strength between geometry and dynamic characteristics. When the comprehensive evaluation result exceeds the preset standard, it indicates that there is a stable intrinsic correlation between trajectory morphology and motion state.

[0023] For trajectory representations with strong coupling relationships, a dimensionality reduction analysis method is used to extract the joint feature directions that contribute most to the overall change, and a low-dimensional feature space is constructed that can simultaneously preserve geometric structure and dynamic information, thereby reducing data dimensionality while retaining key coupling information. This low-dimensional representation is then restored to a new sequence of trajectory sample points and input into a recurrent neural network with a memory mechanism to capture the matching pattern between the curvature change frequency and the velocity change frequency, as well as the trend correspondence between the distribution of curved segment length and the statistical distribution of acceleration amplitude, thus forming a dynamic pattern vector that can characterize the overall dynamic features of the trajectory.

[0024] After obtaining the dynamic mode vector, the correlation between path complexity and energy consumption characteristics is further analyzed. When the correlation is insufficient, the trajectory control points are gradually adjusted through iterative optimization to make the trajectory gradually approach a state with reasonable energy distribution while satisfying the conditions of position continuity, velocity continuity, and acceleration continuity.

[0025] After ensuring kinematic continuity, the trajectory is mapped onto the 3D workspace model of the industrial robot. Combined with the spatial distribution of obstacles, the collision risk caused by the synchronization of bending changes and acceleration peaks within potential interference areas is assessed. When the risk exceeds a set standard, a path search algorithm is used to reconstruct a sequence of candidate nodes to avoid high-risk areas, resulting in a safer obstacle avoidance path.

[0026] Finally, the optimized path undergoes structural consistency verification to ensure that the self-similar geometric structure and the autocorrelation characteristics of the dynamic sequence are not compromised. Simultaneously, a comparative analysis with the original trajectory is used to assess the consistency of the matching degree between key symmetric structures and motion rhythms. If consistency is insufficient, the joint feature weights in the dimensionality reduction stage are adjusted through an error feedback mechanism.

[0027] S1 includes acquiring trajectory data during the operation of an industrial robot, parsing the trajectory data to obtain position, velocity, and acceleration sequences; inputting the position sequence into a spatial feature extraction convolutional neural network to generate a geometric feature vector, and inputting the velocity and acceleration sequences into a dynamic feature extraction convolutional neural network to generate a dynamic feature vector; performing dimensional mapping and feature concatenation on the geometric feature vector and the dynamic feature vector, and fusing them to obtain an initial representation vector containing complete spatiotemporal information.

[0028] In this embodiment, the operation process trajectory data is first obtained from the trajectory recording channel of the industrial robot controller. The trajectory data includes the timestamp, end-effector pose or position, and running status identifier for each sampling moment. The sampling interval is set to 0.004 seconds (this value is directly determined by reading the controller interpolation cycle, which is a fixed configuration item and serves as a unified time reference for the entire process). Before the data enters the parsing process, a data integrity check is performed. The check includes that the timestamp increases monotonically, the sampling interval deviation does not exceed 0.0002 seconds (this upper limit of deviation is jointly determined by the controller communication jitter specification and on-site data collection statistics, and is increased by 0.00005 seconds after calculating the maximum sampling interval deviation in three consecutive shifts of data collection to form a safety margin), and the running status identifier is within the valid range. If there are missing sampling points, the missing length is limited to no more than 3 sampling cycles (this length is determined by the controller cache refresh cycle and the data link instantaneous packet loss limit; exceeding this length will cause unacceptable abrupt changes in speed and acceleration calculations). Within this limit, linear interpolation is used to fill in the missing points. During interpolation, positional continuity is checked for one sampling point before and after the gap. The continuity criterion is the upper limit of single-step displacement. The upper limit of single-step displacement is the maximum speed configured for the robot multiplied by the sampling interval and then multiplied by 1.2 (the coefficient 1.2 is determined by statistically analyzing the ratio of the maximum single-step displacement of historical normal trajectories to the single-step displacement corresponding to the maximum speed, and confirming on the verification data that normal points will not be deleted).

[0029] After integrity processing, the trajectory data is parsed to obtain position, velocity, and acceleration sequences. The position sequence is based on the end-effector position output by the controller. If the output is an end-effector pose, the pose-to-position calculation is performed first. The calculation uses parameters from the same kinematics module as the controller to ensure that the parsed position and the controller's executed position share the same coordinate reference. The coordinate reference is uniformly adopted as the robot's base coordinate system. If the original data is in the workpiece coordinate system or the tool coordinate system, a coordinate transformation is performed. The transformation parameters are taken from the controller's currently active coordinate system configuration. The transformation order is fixed as tool to flange first, then flange to base coordinate system. The transformation result undergoes a round-trip transformation consistency check, and the position hysteresis does not exceed 0.05 mm (this threshold is jointly determined by the controller's coordinate transformation numerical accuracy and the repeatability of the field calibration, selecting an upper limit of 0.03 mm for calibration repeatability and adding a numerical calculation margin of 0.02 mm).

[0030] The velocity sequence is calculated from the position sequence: for each sampling time, the difference between the position at that time and the position at the previous sampling time in three coordinate directions is taken, and divided by the sampling interval to obtain the velocity component for that time period; to ensure that the velocity of the first sampling point is usable, the velocity of the first point is the weighted average of the velocities of the previous three sampling points, with weights set to 0.6, 0.3, and 0.1 according to the temporal proximity (this weight is determined by the noise amplification characteristics of discrete difference at the boundary and the requirement for trajectory smoothness; 0.6 is used to improve the contribution of the nearest interval, and 0.1 is used to suppress the error propagation of the far interval). The acceleration sequence is calculated from the velocity sequence: for each sampling time, the difference between the velocity at that time and the velocity at the previous sampling time in three coordinate directions is taken, and divided by the sampling interval to obtain the acceleration component for that time period; the first acceleration point adopts the same boundary processing strategy as the first velocity point. To reduce sensing noise and high-frequency jitter introduced by interpolation, moving average smoothing was performed on the position, velocity, and acceleration sequences. The window length was set to 5 sampling points (this window length was determined by the requirements for fidelity preservation in the trajectory curvature sensitive area and the requirement for suppressing velocity abrupt changes; among the candidate sets of lengths 3, 5, and 7, a length of 5 simultaneously satisfied the constraints of curvature peak position offset not exceeding 1 sampling point and velocity peak suppression amplitude not less than 30%) on the validation data). The smoothing process generated the output by summing each sample value within the window and dividing by the window length. The window was padded with mirror expansion at both ends of the sequence, with a mirror expansion length of 2 sampling points (this length was determined to be consistent with the window length to avoid mean bias at the boundaries).

[0031] After constructing the three types of sequences, the dual-channel feature extraction is performed. The spatial feature extraction convolutional neural network takes the position sequence as input, with the input tensor arranged in chronological order. It has three channels, each corresponding to one of the three coordinate directions. The network structure uses four layers of one-dimensional convolutions stacked together, with a kernel length of 5, a stride of 1, and symmetrical padding of 2 (the kernel length of 5 and the padding width of 2 are determined by the constraint of maintaining consistency in time length before and after convolution, while the stride of 1 is determined to avoid loss of geometric details due to downsampling). The calculation process for each convolution layer is as follows: a continuous window of length 5 slides along the time axis. The input value at each position within the window is multiplied by the corresponding convolution weight, and the multiplication results are accumulated and a bias is added to obtain the response value at the center of the window. This accumulation process is repeated for all sliding positions to form the output feature sequence of that layer. The number of output channels for the convolutional layers is set to 32, 64, 128, and 128 respectively (this channel configuration is determined by both model latency constraints and feature representation capacity; the latency constraint is set to an inference budget of 2 milliseconds within the on-site control cycle, and the number of channels is selected based on the principle of minimizing validation data error while meeting the budget). Each convolutional layer is followed by batch normalization and nonlinear activation. Batch normalization calculates the mean and variance of the output of the same channel on a mini-batch of samples, and uses this mean and variance to perform linear scaling and translation on the output, stabilizing the channel output distribution. The mini-batch size is set to 64 (this value is determined by the upper limit of GPU memory usage and the stability of mean and variance estimation; 64 is chosen from 32, 64, and 128 to balance stability and resource consumption). Nonlinear activation uses linear rectification, and the calculation process involves subthreshold truncation of each output value, setting values ​​less than 0 to 0 (the threshold of 0 is determined by the definition of linear rectification). After four convolutional layers, global average pooling is performed. Global average pooling accumulates all sampling points along the time axis for each channel and divides by the time length to obtain a summation value for each channel, thus forming a fixed-length geometric feature vector. The length of the geometric feature vector is set to 128 (this length is selected from 64, 128, and 256, and the selection criteria are that the accuracy improvement of trajectory geometric classification on the validation data is not less than 0.8 and the inference latency increment does not exceed 0.2 milliseconds. Finally, 128 was selected to satisfy both constraints).

[0032] The convolutional neural network for dynamic feature extraction takes velocity and acceleration sequences as inputs, with 6 input channels corresponding to the three coordinate components of velocity and acceleration, respectively. This network also employs four layers of one-dimensional convolutions, with a kernel length of 5 and a stride of 1. The padding method is consistent with the spatial channels to ensure temporal alignment. Convolution calculations, batch normalization calculations, linear rectification calculations, and global average pooling calculations are consistent with the spatial channels. Only the input and output channels are configured to enhance high-frequency dynamic features. The output channels are set to 64, 128, 128, and 128 respectively (this configuration is determined by the dynamic feature's requirement to identify velocity and acceleration peaks; within the same inference budget, the first two layers are prioritized to enhance the ability to capture local temporal patterns). The dynamic feature vector length is set to 128 (determined in the same way as the geometric feature vector, requiring both types of vectors to have the same length for subsequent mapping and concatenation).

[0033] To ensure the fusion of the two types of features at the same scale, dimensionality mapping is performed first. Dimensional mapping employs a two-layer fully connected structure to process the geometric shape feature vector and the dynamic feature vector separately. The fully connected computation process is as follows: each component of the input vector is multiplied by its corresponding weight, and the results are accumulated to form the output component. Then, a bias is added to complete the linear transformation from the input space to the target space. Nonlinearity is then introduced through linear rectification. The target dimension after mapping is set to 128 (this dimension is consistent with the length of the aforementioned feature vectors to avoid the inference burden caused by unnecessary dimensionality expansion). The fully connected weights are obtained through offline training. The training data consists of stratified sampling of on-site trajectory data according to the process scenario and the proximity of obstacles. The stratified sampling ratio is set to 1:1:1 (this ratio is determined to avoid the bias risk caused by a few scenarios dominating training; after sampling, the number of samples in each layer is truncated to ensure balance). Training employs a weighted combination of cross-entropy loss and mean squared error loss, with weights set to 0.7 and 0.3 respectively (these weights were determined by evaluating the stability of geometric classification and dynamic regression error on the validation data; 0.7 prioritizes geometric shape recognition, while 0.3 preserves dynamic variation constraints). The learning rate was set to 0.001 (this value was determined by evaluating the rate of loss decrease and oscillation amplitude over five consecutive training epochs; 0.001 achieves the minimum oscillation between the rate of decrease and stability), and the weight decay coefficient was set to 0.0001 (this coefficient was determined by evaluating the overfitting metric; the overfitting metric is set to a difference between the training error and the validation error not exceeding 0.02). Training terminated when the validation error showed no decrease for 10 consecutive epochs (this number of epochs was determined by statistically analyzing the short-term fluctuation cycle of the validation error; the upper limit of the statistical fluctuation cycle was 6 epochs, and 10 epochs were used to exclude occasional fluctuations).

[0034] After mapping, feature concatenation is performed. The concatenation process involves joining two vectors of length 128 end-to-end in a fixed order to form a vector of length 256. The order is fixed: first the geometric feature vector, then the dynamic feature vector (this order is determined by the input specifications of subsequent modules; a fixed order avoids semantic drift in the fused vector). To ensure the fused vector contains complete spatiotemporal information, standardization is performed on both vectors before concatenation. The standardization process involves subtracting the mean of the training data from each component and then dividing by the standard deviation of the training data. The mean and standard deviation of the training data are calculated once during the training phase and saved as model parameters. The lower limit of the standard deviation is set to 0.001 (this lower limit is determined to prevent numerical amplification due to excessively small standard deviations; 0.001 is taken from the floating-point numerical stability constraint and the safety boundary of the minimum standard deviation of the training data, 0.002). The standardized concatenated vector is output as the initial representation vector.

[0035] S2 includes analyzing the initial representation vector to obtain geometric and dynamic feature components, extracting local curvature change sequences and acceleration peak sequences from the geometric and dynamic feature components, calculating the synchronization degree values ​​of the local curvature change sequences and acceleration peak sequences, determining the spatial overlap density of curvature extreme points and velocity abrupt change points based on the synchronization degree values, calculating the directional consistency index based on the spatial overlap density, and determining that the geometric shape and dynamic characteristics have a strong coupling relationship if the comprehensive coupling coefficient generated by the synchronization degree values, spatial overlap density, and directional consistency index meets the conditions.

[0036] In this implementation, the initial representation vector is formed in a fixed splicing order during the generation stage, namely, geometric feature components first, followed by dynamic feature components. The segmentation position is taken as the midpoint of the initial representation vector (the midpoint position is determined by the dimension setting during the feature fusion stage; the dimensions of the geometric feature components and the dynamic feature components are uniformly set to the same value during model finalization to avoid ambiguity caused by inconsistent dimensions after segmentation). At the start of the parsing process, an integrity check is performed on the initial representation vector corresponding to each sampling moment of each trajectory. The check includes that the vector length is consistent with the length recorded in the deployment parameter table, the lengths of the first half and the second half are consistent, and the vector components are finite values ​​and there are no non-numerical terms (the length is directly given by the parameter table generated during model export; the finite value check is used to exclude abnormal components caused by communication jitter and overflow. The check rule adopts component-by-component judgment, and when an abnormal component is found, the sampling moment is marked as invalid to avoid the propagation of abnormal values ​​to the synchronous judgment).

[0037] After segmentation, the calculation process for extracting the local curvature change sequence from the geometric feature components is executed point-by-point. First, the spatial position corresponding to each sampling point is located. This spatial position is taken from the coordinate value aligned with the sampling time in the position sequence obtained from trajectory data parsing (alignment is based on consistency with the timestamp index, which is determined by the sampling interval and the start time and unified in the preprocessing stage). Then, a local three-point neighborhood is constructed for each sampling point. The neighborhood consists of the previous sampling point, the current sampling point, and the next sampling point (the neighborhood span is one sampling point before and after the current sampling point. This span is determined by the sampling interval and the minimum resolvable scale of trajectory geometric change. The minimum resolvable scale is the maximum displacement of adjacent sampling points at the robot's maximum speed, ensuring that local directional changes can be identified within the neighborhood). For this three-point neighborhood, the first segment of the direction vector is calculated, which is the displacement direction from the previous sampling point to the current sampling point; then the second segment of the direction vector is calculated, which is the displacement direction from the current sampling point to the next sampling point; finally, the angle between the two segment of the direction vector is calculated. The larger the angle, the sharper the turn. To transform the angle change into a numerical representation of the local curvature, the angle and the average of two displacement lengths within the neighborhood are used for normalization. Specifically, the average displacement length is used as a scale factor to scale the turning intensity corresponding to the angle, making the same angle comparable at different step sizes (the scale factor is the average of two displacement lengths, used to suppress anomalous amplification of unilateral displacement; the displacement length is obtained from the difference in three-dimensional coordinates and obtained by taking the square root of the sum of squares, calculated by accumulating each item along the coordinate direction). After obtaining the curvature value for each sampling point, they are arranged in chronological order to form a curvature sequence. The local curvature change sequence is obtained by differencing adjacent samples in the curvature sequence. The differencing process involves subtracting the curvature value of the previous sampling point from the current sampling point's curvature value to obtain the curvature change amplitude of that sampling point, and then arranging them in chronological order to form a local curvature change sequence (the difference span is one sampling point, which is determined by the time resolution requirement of the event synchronization judgment, and the time resolution requirement is the sampling interval itself). To suppress the interference of single-point noise on extreme value detection, median filtering is performed on the local curvature change sequence. The filtering window length is set to 3 sampling points (the window length of 3 is determined by the statistical results that impulse noise usually appears in the form of single or double points, while satisfying the constraint that the extreme value position offset does not exceed 1 sampling point; this constraint is verified on the labeled data, and the verification standard is that the maximum deviation between the curvature event position and the manually labeled event position does not exceed 1 sampling point).

[0038] The calculation process for extracting the acceleration peak sequence from the dynamic feature components is also performed point-by-point. First, the acceleration vector corresponding to each sampling point is located. This acceleration vector is taken from the three coordinate direction components aligned with the sampling time in the acceleration sequence obtained from trajectory data parsing (alignment is based on consistent timestamp indices). Then, the acceleration amplitude of each sampling point is calculated. The acceleration amplitude is obtained by summing the squares of the three direction components and taking the square root. The calculation is completed by accumulating each component in each direction. Peak detection is performed on the acceleration amplitude sequence. Peak detection uses a sliding window with a length of 5 sampling points. The window is centered on the current sampling point and consists of the two sampling points before and the two sampling points after it (the window length of 5 is determined by the maximum number of lag sampling points from the servo system command to the feedback. The maximum number of lag sampling points is obtained by statistically analyzing the alignment error between the command and feedback of multiple acceleration and deceleration trajectories on site. The maximum value of the statistical result is taken as 2 sampling points, and then the center point is incorporated to form a length of 5). Within each window, if the acceleration amplitude at the center point is strictly greater than the acceleration amplitudes of the other four sampling points within the window, the center point is marked as a peak point and its amplitude is output; otherwise, 0 is output. These peak values ​​are arranged in chronological order to form an acceleration peak sequence. To avoid excessively dense peaks causing event count bias in synchronization determination, a minimum interval constraint is introduced for the peak points. The minimum interval is set to 3 sampling points (the minimum interval of 3 is jointly determined by the controller interpolation period and the shortest significant acceleration change interval achievable by the mechanical system; the minimum interval between adjacent true peaks on the labeled data is 3 sampling points, and this value is used to eliminate false peaks caused by noise).

[0039] The synchronization degree value is used to characterize the alignment strength between the local curvature change sequence and the acceleration peak sequence on the time axis. Its calculation is based on event matching. First, curvature events are detected from the local curvature change sequence. The detection rule adopts a sliding window with a length of 5 sampling points. The window is centered on the current sampling point and consists of the first 2 sampling points and the last 2 sampling points. When the curvature change amplitude of the center point is strictly greater than that of the other 4 sampling points in the window, the center point is marked as a positive curvature event. When the curvature change amplitude of the center point is strictly less than that of the other 4 sampling points in the window, the center point is marked as a negative curvature event (the window length of 5 is consistent with the acceleration peak detection window to ensure that the two types of events have consistent time resolution; both positive and negative events are retained to cover both curvature intensification and curvature deterioration changes). Subsequently, based on each curvature event, the time search range is defined as two sampling points before and after the event (the two sampling points before and after the search range are determined by the maximum value of the alignment error between the instruction and the feedback, and a safety margin of one sampling point is added to the maximum value to form the final range; the safety margin is used to cover the offset of one sampling point caused by occasional communication jitter). Within the time search range, the peak acceleration point is searched. If a peak point exists, the peak point with the smallest time offset is selected as the matching point, and a matching score is given based on the time offset: 1.0 for a time offset of 0 sampling points (full score is given for perfect alignment); 0.7 for a time offset of 1 sampling point (the value is determined by the alignment error distribution of strongly coupled samples in the labeled data; samples with an error of 1 sampling point account for a significant proportion, and the score must be higher than 0.5 to maintain stable identification of strong coupling); 0.4 for a time offset of 2 sampling points (the value is determined by the fact that strongly coupled samples with an error of 2 sampling points still exist in the labeled data but their proportion has decreased; the score must be lower than 0.5 to increase the distinguishability from those with errors of 0 or 1). If no peak point exists within the time search range, the matching score for the curvature event is 0. The synchronization score is obtained by summing the matching scores of all curvature events and dividing by the number of curvature events (the summation process is to add up the events one by one, and the division uses the number of events as the denominator to form the average; when the number of curvature events is 0, the synchronization score is directly taken as 0 to avoid the denominator being 0).

[0040] Spatial overlap density is used to characterize the spatial intensity of the correspondence between curvature extrema and velocity abrupt change points. Velocity abrupt change points are calculated from the velocity sequence and are constrained by the degree of synchronization. First, the velocity amplitude is calculated from the velocity sequence obtained by parsing the trajectory data. The velocity amplitude is obtained by taking the square root of the sum of the squares of the three directional velocity components, and the calculation is performed by accumulating each component in each direction. Then, the velocity variation amplitude sequence is calculated. The velocity variation amplitude sequence consists of the absolute values ​​of the differences between the velocity amplitudes of adjacent sampling points. Specifically, for each sampling point, the current velocity amplitude is subtracted from the velocity amplitude of the previous sampling point, and the absolute value is taken, arranged in chronological order to form the velocity variation amplitude sequence. Speed ​​abrupt change points are detected using threshold discrimination. The threshold is set as the maximum value of the speed change amplitude sequence within the normal operating segment of the labeled data plus a fixed safety margin (the normal operating segment is defined as a segment without emergency stop or acceleration commands, and the segment selection is based on the controller's operating status identifier; the maximum value is obtained by traversing all normal segments point by point; the safety margin is 0.1 times this maximum value, which is estimated from the worst-case scenario after superimposing the upper bound of speed measurement noise and the upper bound of interpolation error, and this multiple is used to avoid noise-triggered abrupt changes). When the speed change amplitude of a sampling point is strictly greater than this threshold, the sampling point is marked as a speed abrupt change point. Subsequently, spatial coincidence determination is performed, and the spatial matching tolerance is determined by a combination of the basic tolerance and the synchronization scaling factor. The basic tolerance is the sum of the robot's repeatability accuracy and the sampling spatial resolution (the repeatability accuracy is the maximum repeatability error given in the equipment calibration report; the sampling spatial resolution is the maximum single-step displacement obtained by multiplying the maximum speed by the sampling interval, the maximum speed is the speed limit value in the controller's safety parameters, and the sampling interval is the trajectory data sampling period). The synchronization scaling factor is set in three levels based on the synchronization level: 0.8 is used when the synchronization level is not lower than 0.8 (the 0.8 threshold is determined by the lower bound of the synchronization level distribution of strongly coupled samples in the labeled data; the scaling factor of 0.8 is used to tighten the tolerance and reduce false matches); 1.0 is used when the synchronization level is between 0.5 and 0.8 (the 0.5 threshold is determined by the upper bound of the synchronization level distribution of weakly coupled samples; 1.0 is used to maintain the basic tolerance); and 1.2 is used when the synchronization level is lower than 0.5 (1.2 is used to loosen the tolerance to cover spatial offsets when time alignment is poor; the value is determined by constraining the number of false matches on the labeled data to not exceed 0.2 of the number of matches and maximizing the number of effective matches). The spatial matching tolerance is the basic tolerance multiplied by the synchronization scaling factor. For each curvature event sampling point, a candidate velocity mutation point set is constructed. The candidate set is limited to the time range of two sampling points before and after the curvature event (the time range is consistent with the synchronization matching to ensure that the events are compared under the same response lag scale). For each velocity mutation point in the candidate set, calculate the straight-line distance between its spatial position and the spatial position of the curvature event. The straight-line distance is obtained by taking the square root of the sum of the squares of the differences in the three coordinate directions. The calculation is completed by accumulating the values ​​in each direction.If a candidate point exists whose straight-line distance is less than or equal to the spatial matching tolerance, the candidate point with the smallest straight-line distance is selected as the coincident matching point, and the curvature event is recorded as 1 coincidence. If no candidate point satisfies the tolerance, it is recorded as 0 coincidences. To avoid the same velocity mutation point being repeatedly occupied by multiple curvature events, a one-to-one matching rule is adopted. Once a velocity mutation point is selected, it is removed from the subsequent candidate set. The spatial coincidence density is obtained by summing the coincidence counts of all curvature events and dividing by the number of curvature events (the spatial coincidence density is 0 when the number of curvature events is 0).

[0041] The directional consistency index is used to characterize the degree of consistency between the direction of geometric change and the direction of dynamic change. Its calculation uses matched pairs of coincident points as input. For each pair of coincident points, the direction of geometric change is first calculated. This direction is based on the tangential change of the trajectory: taking the preceding and following sampling points of the coincident points to form two adjacent trajectory segments, the motion directions of the preceding and following segments are calculated separately. The motion direction is obtained by normalizing the displacement vectors of adjacent points. The normalization process involves dividing each directional component of the displacement vector by the displacement length, which is obtained by taking the square root of the sum of squares. Subsequently, the difference vector between the preceding and following motion directions is used as the tangential change direction, obtained by subtracting the corresponding components. Next, the direction of dynamic change is calculated. This direction is taken as the direction of the difference vector between the acceleration vector of the coincident point and the acceleration vector of the preceding sampling point, obtained by subtracting the corresponding components. Then, normalization is performed according to the length of the difference vector, using the same method as for the displacement vector. The angle between the two directions is then calculated. This angle is obtained by the dot product of the direction vectors and their respective lengths. The dot product is obtained by multiplying the corresponding components and summing them. The angle is mapped to a consistency score using fixed segments: an angle no greater than 15 degrees scores 1.0 (15 degrees is determined by the upper bound of the angle distribution of strongly coupled samples, ensuring that strongly coupled samples fall within the full-score segment); an angle greater than 15 degrees but no greater than 45 degrees scores 0.7 (45 degrees is determined by the boundary between the angle distributions of strongly coupled and weakly coupled samples, used to form a stable distinction); an angle greater than 45 degrees but no greater than 90 degrees scores 0.4 (90 degrees corresponds to orthogonality, 0.4 is used to retain some weak consistency contribution and differentiate it from strong consistency); an angle greater than 90 degrees scores 0 (greater than 90 degrees corresponds to opposite directions, 0 is used to explicitly negate consistency). The consistency scores of all overlapping point pairs are summed and divided by the number of overlapping point pairs to obtain the directional consistency index (the directional consistency index is 0 when the number of overlapping point pairs is 0).

[0042] The comprehensive coupling coefficient is obtained by weighting three factors: synchronization degree, spatial overlap density, and directional consistency index. The weights are taken from the three coefficients of a linear combination of the three factors, and the sum of the coefficients is 1.0 (the sum of the coefficients of 1.0 is used to keep the numerical range of the comprehensive coupling coefficient stable and facilitate the standardization of the threshold). The three weights are determined through the parameter selection process of the labeled data: First, a candidate weight set is set, with candidate weights ranging from 0.1 to 0.8 and a step size of 0.1 (the step size of 0.1 is determined by the constraints of the number of candidate combinations and the upper limit of the verification calculation, ensuring that the full combination evaluation is completed within a single model review cycle); for each weight combination, the comprehensive coupling coefficient is calculated for each trajectory and compared with the strong coupling conclusion of the label, and the detection rate and false positive rate are statistically analyzed; under the constraint that the false positive rate is not higher than 0.1 and the detection rate is not lower than 0.9, the weight combination with the highest accuracy is selected (the upper limit of the false positive rate of 0.1 and the lower limit of the detection rate of 0.9 are determined by the safety constraints of the obstacle avoidance scenario, and the detection rate is improved under the premise of controlled false positives to avoid missed detections). The threshold for determining the overall coupling coefficient is also determined using labeled data: A scan threshold with a step size of 0.01 is used within the range of 0.00 to 1.00 (the step size of 0.01 is determined by the combined constraints of the overall coupling coefficient resolution requirement and scanning cost). For each threshold, the detection rate and false positive rate are statistically analyzed, and the accuracy is calculated. The threshold with the highest accuracy among those satisfying a detection rate of not less than 0.9 and a false positive rate of not more than 0.1 is selected as the final threshold. When the overall coupling coefficient is greater than or equal to this threshold, a strong coupling relationship is determined between the geometry and dynamic characteristics.

[0043] S3 includes obtaining the initial representation vector selected by strong coupling judgment, constructing a joint covariance matrix using the initial representation vector; decomposing the joint covariance matrix to extract the geometric dynamic coupling principal direction that retains the maximum variance, and identifying the joint principal component axis in the cross-modal information set; constructing an initial subspace based on the geometric dynamic coupling principal direction and the joint principal component axis, and optimizing the target subspace by minimizing the joint reconstruction error function; mapping the initial representation vector to the target subspace to obtain the dimensionality-reduced joint feature representation.

[0044] In this implementation, initial representation vectors marked as strongly coupled are first extracted from the strong coupling determination results, and then the vectors are grouped into the same trajectory segment according to the time sequence of the same trajectory. The length of the trajectory segment is set to 256 sampling points (this length is determined by the controller interpolation cycle and the upper limit of the on-site inference latency: the interpolation cycle is directly read from the controller configuration, and the upper limit of latency is determined by the time budget that the trajectory planning module must complete the calculation within a single planning cycle; among the candidate lengths of 128, 192, 256, and 320, when the length is 256, a single segment covers the continuous interval of common acceleration, deceleration, and steering combination actions, while maintaining a stable single-segment dimensionality reduction calculation latency below the upper limit on the target hardware). Each sampling point corresponds to one initial representation vector, and the vector length is set to 256 (this length is determined by the output structure of the feature fusion stage, with 128 geometric feature components and 128 dynamic feature components, which are concatenated in a fixed order to obtain 256; the deployment side reads this value through the model parameter table to avoid human configuration differences). Before proceeding with covariance statistics, three checks are performed on each vector sequentially: First, length check, verifying that the number of components is 256 for each vector; second, numerical check, checking for non-numerical terms or infinity for each component; third, amplitude check, checking whether the amplitude exceeds the upper limit of the component amplitude calculated during the training phase. The upper limit of amplitude is the maximum absolute value of the corresponding component in the training data plus a fixed margin of 0.05 (the margin of 0.05 is determined by combining the upper limit of instantaneous deviation introduced by floating-point quantization error and communication jitter, and is chosen to ensure that normal data is not mistakenly removed and abnormal spikes are intercepted). Any vector that fails any check is removed from the sample set. The removal is done by marking it as invalid in situ and skipping subsequent calculations to avoid outliers participating in covariance accumulation and causing a shift in the main direction. The minimum size of the sample set is set to 512 vectors (this size is determined by the stability of covariance estimation: during the model finalization stage, the fluctuation of the principal direction angle under different sample sizes is statistically analyzed. When the size is less than 512, the fluctuation of the principal direction between adjacent batches increases significantly, which leads to instability in the dimension reduction projection. After the size is set to 512, the fluctuation enters the stable range, while satisfying the constraints of on-site buffer capacity and processing latency).

[0045] The construction of the joint covariance matrix follows the sequence of "calculating the mean dimension by dimension, removing the mean from each sample, accumulating the out-of-sample values, and normalizing by sample size." First, the mean vector is calculated dimension by dimension for the sample set: for the first dimension, the values ​​of all samples in that dimension are summed, and then divided by the number of samples to obtain the mean for that dimension; the same summation and division process is repeated for the second to 256th dimensions to form a mean vector of length 256. Then, the mean is removed from each sample vector: for the first dimension of the sample vector, the first dimension of the mean vector is subtracted to obtain the centered component; the second to 256th dimensions are subtracted sequentially to obtain the centered sample vector. Next, the outer product is calculated and accumulated: For the centered sample vector, the first-dimensional component is multiplied by itself to obtain the element in the first row and first column of the matrix; the first-dimensional component is multiplied by the second-dimensional component to obtain the element in the first row and second column of the matrix, and so on, until the first-dimensional component is multiplied by the 256th-dimensional component to obtain the element in the first row and 256th column of the matrix; then the second-dimensional component is multiplied by the first-dimensional component to obtain the element in the second row and first column of the matrix, and so on, until the element in the 256th row and 256th column is obtained, forming a 256-row, 256-column outer product matrix. This outer product matrix is ​​then accumulated element by element into the accumulation matrix, which is initially all zeros. The accumulation process is repeated for each centered sample in the sample set. Finally, scaling is performed on the cumulative matrix: each element of the cumulative matrix is ​​divided by the number of samples minus 1 to obtain the joint covariance matrix (the denominator is the number of samples minus 1 to reduce the systematic bias of underestimating variance in a finite number of samples; this value is determined after statistical verification of variance bias from multiple samplings of the same data during the model finalization stage). To improve the stability of the decomposition, diagonal stabilization is performed on the joint covariance matrix: a stabilization increment of 0.000001 is added to the element in the first row and first column, and the same increment is added to each item from the second row and second column to the 256th row and 256th column (the increment of 0.000001 is jointly determined by the floating-point precision and the risk of decomposition failure when the matrix is ​​close to singular; during the finalization stage, multiple increments are scanned with the constraint of "the main direction angle offset is no more than 1 degree and there is no numerical collapse in the decomposition process", and finally 0.000001 is selected to satisfy both constraints and keep the reconstruction error from deteriorating).

[0046] The joint covariance matrix is ​​decomposed using a symmetric matrix decomposition process to obtain a set of orthogonal directions sorted by variance contribution. The decomposition process employs an iterative solution: each iteration updates the direction set and corresponding variance contribution values, and calculates the maximum change in variance contribution between two adjacent iterations; iteration stops when the maximum change is no higher than 0.0000001 (the threshold of 0.0000001 is determined by the magnitude of the matrix elements and the stability of the dimensionality reduction output; during the modeling phase, this threshold is correlated with the fluctuation of the projected coordinates statistically; when the threshold is 0.0000001, the fluctuation of the projected coordinates under repeated calculations enters the noise-dominated range, and continued iteration does not bring verifiable improvement). The maximum number of iterations is set to 200 (this upper limit is determined by the statistics of the slowest convergence batch; the number of iteration convergence steps is recorded on multiple batches of field trajectories, 200 covers all batches and ensures that the worst-case latency remains under control). After decomposition, the variance contribution values ​​are sorted from largest to smallest, and the corresponding direction vectors are rearranged synchronously. The first direction vector after sorting is defined as the geometric-dynamic coupling principal direction (this direction carries the maximum discrete information in the sample set and satisfies the definition of the principal direction of "maximum variance preservation"). The identification of the joint principal component axis is completed in two layers: the first layer determines the number of principal components, and the second layer selects axes from the principal components that are in the cross-modal information set. The number of principal components is determined by the variance explanation ratio threshold of 0.9 (the threshold of 0.9 is determined by the compression ratio, reconstruction error and subsequent trajectory reconstruction availability constraints: during the modeling stage, the thresholds of 0.85, 0.88, 0.90, 0.92 and 0.95 are evaluated one by one, requiring that the dimensionality reduction dimension does not exceed 32 and the reconstruction error decreases significantly while maintaining the feasibility rate of subsequent trajectory optimization remains stable. The final threshold of 0.9 satisfies the three constraints). The calculation of the variance proportion is performed by dividing the cumulative variance contribution by the total variance contribution: starting from the first variance contribution value after sorting, each value is accumulated to obtain the cumulative value; all 256 variance contribution values ​​are accumulated to obtain the total value; the cumulative value is divided by the total value to obtain the proportion; accumulation stops when the proportion reaches or exceeds 0.9, and the number of directions accumulated at this point is determined as the number of principal components. The screening across modal axes uses the modal energy balance criterion: for each candidate direction vector, the first 128 components are taken as geometric components, and the last 128 components are taken as dynamic components; the geometric energy is obtained by squaring each geometric component and accumulating them, and the dynamic energy is obtained by squaring each dynamic component and accumulating them; the smaller of the two energies is taken as the lower bound energy, and the lower bound energy is divided by the sum of the two energies to obtain the balance score. The equilibrium score threshold was set to 0.35 (the threshold of 0.35 was determined during the modeling phase using labeled data: the equilibrium score distribution of candidate directions for strongly coupled and weakly coupled samples was statistically analyzed, and the threshold that minimized the overlap between the two types of samples and achieved the highest accuracy under the constraints of "false positive rate not higher than 0.1 and detection rate not lower than 0.9" was selected, ultimately set to 0.35). Direction vectors with an equilibrium score not lower than 0.35 were labeled as joint principal component axes and sorted in descending order of variance contribution.To control the size of the subspace, the upper limit for retaining the joint principal component axes is set to 16 (this upper limit is determined by the matrix multiplication size and inference delay budget: projection calculation requires component-wise multiplication and addition of each vector and each axis, and the more axes there are, the greater the number of multiplication and addition; when the upper limit is 16, the cross-modal information is sufficiently retained and the single-segment processing delay is kept within the budget).

[0047] The initial subspace is constructed using the following steps: axis set assembly, redundancy removal, orthogonalization, and dimension pruning. During axis set assembly, the geometrically and dynamically coupled principal directions are first added, followed by the principal component axes corresponding to the aforementioned number of principal components, along with no more than 16 joint principal component axes. Redundancy removal employs a directional proximity criterion: for any two axes, their included angle is calculated. This angle is based on the dot product of the two axes and their respective lengths. The dot product is obtained by multiplying the corresponding components one by one and then summing them. The lengths are obtained by summing the squares of the components and then taking the square root. An angle less than 5 degrees is considered directionally close, and only the axis with the larger variance contribution is retained (the 5-degree threshold is determined by the numerical error and the upper bound of the projection difference: during the finalization stage, statistics show that when the included angle is within 5 degrees, the projection difference between the two axes contributes less to the reconstruction error than the noise fluctuation amplitude, so only one axis needs to be retained). Orthogonalization employs axis-by-axis projection elimination: the current axis is processed sequentially according to the axis set. First, the projection component of the current axis on the space spanned by the processed axes is calculated. The projection component is obtained by "obtaining coefficients through inner product, multiplying the coefficients by the basis axes, and then subtracting each term". After deducting the projection component, the remaining component orthogonal to the processed axis is obtained. Then, the remaining component is normalized to a unit length according to its length. Normalization is completed by dividing each component by its length. After orthogonalization, dimension pruning cuts the axis set into 32 axes to form an initial subspace (the dimension of 32 is jointly determined by three constraints: compression ratio, reconstruction error, and the input size of subsequent temporal modeling: during the shaping stage, 16, 24, 32, and 40 are evaluated one by one, requiring that the length of the vector after dimension reduction does not exceed 32 and the reconstruction error of subsequent trajectory sample points does not exceed the preset upper limit while maintaining the overall time delay budget. The final dimension is 32 to satisfy the constraints and minimize the stable reconstruction error).

[0048] The optimization of the target subspace aims to minimize the joint reconstruction error. The optimization process is executed according to the steps of "projection reconstruction, residual calculation, joint error summarization, subspace update, and termination determination," maintaining a dimension of 32 in each iteration. In the projection reconstruction stage, for each sample vector, its projected coordinates on the 32 basis axes of the current subspace are calculated sequentially: for the first basis axis, the inner product of the sample vector and this basis axis is calculated, and the inner product is multiplied by the corresponding components and then summed; this multiplication and summation process is repeated for the second to the 32nd basis axes to obtain 32 coordinate values. Then, the reconstruction vector is reconstructed based on the coordinate values ​​and the basis axes: for the first basis axis, each component of this basis axis is multiplied by the first coordinate value to obtain the contribution vector of that axis; for the second to the 32nd basis axes, the corresponding contribution vectors are obtained; the 32 contribution vectors are summed by their component positions to form the reconstruction vector. In the residual calculation stage, the original sample vector is subtracted from the reconstruction vector component by component to obtain the residual vector, and the first 128 dimensions of the residual vector are defined as the geometric residual, and the last 128 dimensions are defined as the dynamic residual. In the joint error summarization stage, the mean square error of the geometric residuals and the mean square error of the dynamic residuals are calculated separately: the geometric residuals are squared component by component and summed, then divided by 128 to obtain the geometric mean square error; the dynamic residuals are squared component by component and summed, then divided by 128 to obtain the dynamic mean square error. To avoid a large variance in one mode dominating the optimization, the two types of mean square errors are weighted according to their inverse variance ratio and then summed in a weighted manner: first, the total variances of the geometric and dynamic components are statistically analyzed in the sample set. The variance statistics are performed by "removing the mean from each dimension, squaring and summing, then dividing by the sample size minus 1", and then averaging the variances of each dimension to obtain the total variance; then, the reciprocals of the two total variances are taken as the original weights and normalized so that the sum of the two weights is 1; then, the normalized weights are multiplied by the two types of mean square errors and summed to obtain the joint reconstruction error. In the subspace update phase, a residual covariance matrix is ​​constructed to capture the main directions that are not explained in the current subspace: the mean of all residual vectors is recalculated, the mean is removed, the outer accumulation is performed, the scale is normalized, and the diagonal is stabilized according to the aforementioned covariance construction process to obtain the residual covariance matrix; symmetric decomposition is performed on the residual covariance matrix, and the direction with the largest variance contribution is taken as the update direction; the update direction is incorporated into the current 32 basis axes, redundancy removal and orthogonalization are performed, and the first 32 axes are truncated to form the updated subspace. The termination decision adopts a dual constraint of the joint reconstruction error reduction threshold and the maximum number of iterations: termination occurs when the joint reconstruction error reduction of two adjacent iterations is not higher than 0.00001 (the threshold of 0.00001 is determined by the upper bound of error fluctuation caused by quantization noise in the field data: during the modeling stage, the natural fluctuation amplitude of the joint error is repeatedly calculated and statistically analyzed for the same trajectory. If the fluctuation amplitude stably falls within 0.00001, further optimization will not produce a verifiable reduction); the maximum number of iterations is set to 50 (this upper limit is determined by the iteration delay budget of the worst batch. When the upper limit is 50, the processing of a single segment still does not exceed the budget of the planning cycle).

[0049] After the target subspace is determined, the dimensionality reduction mapping is performed according to the "projected coordinate output and scale shaping" process. For each initial representation vector, the inner product is calculated with each of the 32 basis axes of the target subspace, resulting in 32 projected coordinates. The coordinate order is consistent with the basis axis order and remains fixed. To ensure consistent input scale for subsequent modules, the 32 coordinates are standardized. Standardization uses the coordinate mean and standard deviation obtained during the training phase. The coordinate mean is obtained by summing the projected coordinates of the training samples one by one and dividing by the number of samples. The coordinate standard deviation is obtained by square-scumming the squares of the projected coordinates after removing the mean, dividing by the number of samples minus 1, and then taking the square root. The lower limit of the standard deviation is set to 0.001 (the lower limit of 0.001 is determined by numerical stability constraints: when the standard deviation is less than 0.001, standardization will amplify small noise and form spikes in the input during subsequent time series modeling; when the lower limit is 0.001, the spikes disappear and the coordinate distribution remains separable). The standardized 32-dimensional coordinate sequence serves as the joint feature representation after dimensionality reduction.

[0050] S4 includes obtaining the dimensionality-reduced joint feature representation, reconstructing the industrial robot trajectory sample point sequence through inverse projection; calculating the frequency matching degree of curvature and velocity zero-crossing frequency, and the distribution matching trend of bending segment length and acceleration amplitude based on the industrial robot trajectory sample point sequence; combining the frequency matching degree and distribution matching trend to generate a temporal feature sequence, and inputting the temporal feature sequence into a long short-term memory network; obtaining the hidden layer state vector output by the long short-term memory network, and mapping the hidden layer state vector to obtain the trajectory dynamics mode vector.

[0051] In this embodiment, the joint feature representations after dimensionality reduction are arranged in the order of sampling time to form a sequence. Each sampling time in the sequence corresponds to a set of dimensionality-reduced coordinates. The length of the dimensionality-reduced coordinates is set to 32 (the value of 32 is determined by information retention and computation delay: the reconstruction error and single-segment processing delay are evaluated simultaneously among the candidate lengths of 16, 24, 32, and 40. When the length is 32, the reduction of the reconstruction error reaches a stable range and the single-segment processing delay does not exceed the budget of the planning period). The arrangement order of each set of dimensionality-reduced coordinates is consistent with the order of the target subspace basis axes. The target subspace basis axes are obtained from the previous dimensionality reduction process and are saved along with the model parameters. After obtaining the dimensionality-reduced coordinates, scale recovery is performed first. The mean and standard deviation used for scale recovery are taken from the statistical parameters of the dimensionality reduction stage. The mean is obtained by summing the modeling data coordinates one by one under the same base axis coordinates and then dividing by the number of samples. The standard deviation is obtained by square-squaring each coordinate after removing the mean, summing the squares, dividing by the number of samples minus 1, and then taking the square root. The lower limit of the standard deviation is set to 0.001 (the value of 0.001 is determined by numerical stability: when the standard deviation is less than 0.001, scale recovery will amplify quantization noise and form peak inputs in subsequent time series features. After the lower limit is set to 0.001, the peaks disappear and the coordinate distribution remains separable). Scale recovery is performed item by item for each coordinate: each coordinate is first multiplied by the corresponding standard deviation and then added to the corresponding mean to obtain the recovered coordinates. The recovered coordinates are used for subsequent inverse projection operations.

[0052] The inverse projection reconstruction forms a restored representation vector based on the target subspace basis axes. The length of the restored representation vector is set to 256 (the value 256 is determined by the feature fusion output structure: 128 geometric feature components and 128 dynamic feature components are concatenated in a fixed order to form 256; the deployment side verifies this length using the model parameter table to ensure consistency). The inverse projection calculation is completed by accumulating along each basis axis: for the first basis axis, the first value of the restored coordinate is taken and multiplied by each of the 256 components of that basis axis to obtain the contribution vector of the first basis axis; for the second to the 32nd basis axes, the component-by-component multiplication is repeated to obtain the corresponding contribution vectors; then, the 32 contribution vectors are added one by one at the same component positions to obtain the centered reconstruction vector. Subsequently, the mean vector used for decentralization in the dimensionality reduction stage is added back to the reconstructed vector component by component to obtain the recovered representation vector. The mean vector is obtained by accumulating the initial representation vectors of strongly coupled samples dimension by dimension and then dividing by the number of samples (the mean vector is determined based on the statistical consistency of covariance to ensure that the recovered representation returns to the same central position as the covariance construction, avoiding the systematic reconstruction bias caused by center drift).

[0053] The sequence of industrial robot trajectory sample points is obtained by decoding the sequence of recovered representation vectors. The decoding process is performed point by point according to the sampling time while maintaining consistency of the time index. First, the recovered representation vector at each sampling time is mapped to a three-dimensional position sample point. The mapping adopts a linear decoding structure with fixed weights: for the first coordinate component of the position, the 256 components of the recovered representation vector are multiplied by their corresponding weights one by one and then accumulated, and then the bias is added to obtain the output. The same multiplication and addition process is repeated for the second and third coordinate components of the position to obtain the output, thus obtaining the three-dimensional position sample point at that sampling time. The decoding weights and biases are determined by supervised learning during the modeling stage. The learning objective is to minimize the deviation between the decoded position sample points and the known trajectory positions. The learning data is taken from strongly coupled trajectory segments. The learning process minimizes the mean square loss of the decoding error. The training termination condition is that the verification error does not decrease for 10 consecutive rounds (the value of 10 is determined by the statistical analysis of the short-term fluctuation period of the verification error, with an upper limit of 6 rounds for the fluctuation period, and 10 rounds are used to exclude occasional fluctuations). After decoding, a sequence of position sample points arranged in chronological order is obtained; the sampling interval is set to 0.004 seconds (the value 0.004 is determined by reading the controller interpolation cycle, which is a controller configuration item, and the reading result serves as a unified time reference for the entire process). Velocity and acceleration are calculated based on the position sample point sequence. The velocity vector is calculated by subtracting the position of the previous sample point from the position of the second sample point to obtain the displacement vector. Then, each component of the displacement vector is divided by 0.004 seconds to obtain the velocity vector. The velocity vector of the first sample point is a weighted combination of the first and second segment velocity vectors, with weights set to 0.6 and 0.4 (the values ​​0.6 and 0.4 are determined by boundary differential noise suppression requirements to ensure that the velocity at the first point is not excessively affected by single-segment errors and remains continuous with subsequent velocities). The acceleration vector is calculated using velocity vector differentiation: from the second to the last sampling point, the current velocity vector is subtracted from the previous velocity vector to obtain the velocity change vector. Then, each component of the velocity change vector is divided by 0.004 seconds to obtain the acceleration vector. The acceleration vector at the first sampling point uses the same boundary processing strategy as the initial velocity point. The velocity amplitude is obtained by summing the squares of the three components of the velocity vector at each sampling point and taking the square root. Similarly, the acceleration amplitude is obtained by summing the squares of the three components of the acceleration vector at each sampling point and taking the square root. This forms a sequence of trajectory sample points that includes position, velocity amplitude, and acceleration amplitude and is time-aligned.

[0054] The curvature zero-crossing frequency and velocity zero-crossing frequency are calculated on the same sliding time window. The time window length is set to 64 sampling points, and the sliding step size is set to 16 sampling points (the values ​​64 and 16 are determined by both time resolution and statistical stability: after evaluating the zero-crossing count variance and response delay to sudden changes, the candidate lengths of 32, 64, and 96 show that the count variance decreases significantly and the response delay meets the planning period when the length is 64; the step size of 16 allows adjacent time windows to share 48 sampling points, ensuring a smooth frequency sequence and a computational cost lower than the processing budget). Curvature calculation is performed point by point: for the second to the second-to-last sampling point, a three-point neighborhood is formed by the previous sampling point, the current sampling point, and the next sampling point, and the displacement vectors of the previous and next segments are calculated respectively; each component of the displacement vector is obtained by subtracting the corresponding coordinates, and the displacement length is obtained by summing the squares of each component and taking the square root; each component of the displacement vector is divided by the displacement length to obtain the unit direction vector. The angle between the two unit direction vectors is then calculated. This angle is derived from the dot product of the two vectors and their respective lengths. The dot product is obtained by multiplying the corresponding components and summing them. The curvature amplitude is obtained by normalizing the angle and the average neighborhood step size. The average neighborhood step size is the arithmetic mean of the lengths of the preceding and following displacement segments (this average is used to suppress abnormal curvature amplification by unilateral step size). The curvature sign is determined using a consistent steering discrimination rule: using a preset vertical direction in the robot's base coordinate system as the reference direction, the steering orientation formed by the three-point neighborhood is calculated. When the steering orientation is in the same direction as the reference direction, the curvature sign is positive; when the steering orientation is in the opposite direction, the curvature sign is negative. The steering orientation is determined by the normal direction of the steering plane formed by the preceding and following directions. This normal direction is determined by the steering relationship between the two directions in three-dimensional space. Using the same reference direction ensures consistent signs for the same trajectory segment. Curvature zero-crossing points are detected by sign change: traversing the curvature sign along the time axis, if two adjacent sampling points have opposite curvature signs and both have curvature amplitudes higher than the lower limit of curvature noise, then a zero-crossing is determined between the two points. The lower limit of curvature noise is determined during the modeling phase: a trajectory segment labeled as linear motion is selected, requiring that the acceleration amplitude of this segment be continuously lower than the steady-state threshold for a duration of no less than 32 sampling points. The steady-state threshold is the upper limit of servo steady-state acceleration jitter multiplied by 1.2 (the multiplier 1.2 is estimated and determined by the worst-case scenario of the jitter upper limit and the sampling jitter superimposed); the maximum curvature amplitude is counted within the straight segment that meets the conditions, and then multiplied by 1.1 to obtain the lower limit of curvature noise (the multiplier 1.1 is used to cover batch-to-batch differences and avoid occasional noise in straight segments triggering zero-crossing). The curvature zero-crossing frequency is calculated within each time window: the number of zero-crossings is counted within the sampling point range covered by the window, and the frequency of the window is obtained by dividing the number of zero-crossings by the time window coverage duration. The time window coverage duration is determined by the time length corresponding to the 64 sampling points, i.e., it is obtained by accumulating 63 sampling intervals. The curvature zero-crossing frequency sequence is calculated sequentially for all time windows.

[0055] The velocity zero-crossing frequency is calculated based on the velocity deviation sequence, which is formed by averaging within a time window. For each time window, the velocity amplitudes of the 64 sampling points within the window are first accumulated and then divided by 64 to obtain the velocity mean. Then, the velocity deviation is obtained by subtracting the velocity mean from the velocity amplitude of each sampling point within the window. A positive velocity deviation indicates that it is higher than the mean, and a negative deviation indicates that it is lower than the mean. Velocity zero-crossing points are detected by changes in the sign of the deviation: traversing the velocity deviation sign along the time axis, if two adjacent sampling points have opposite deviation signs, and the deviation amplitudes of both points are higher than the lower limit of velocity noise, then a zero-crossing is determined to have occurred between the two points. The lower limit of velocity noise is determined during the modeling phase: A trajectory segment labeled as uniformly moving is selected, requiring that the acceleration amplitude of this segment be continuously lower than the steady-state threshold for a duration of at least 32 sampling points. The steady-state threshold is the upper bound of the servo steady-state acceleration jitter multiplied by 1.2 (the multiplier 1.2 is determined based on the same criteria as before). Within the uniformly moving segment that meets the conditions, the maximum value of the velocity deviation amplitude is statistically analyzed, and then multiplied by 1.1 to obtain the lower limit of velocity noise (the multiplier 1.1 is used to cover inter-batch differences and suppress noise-triggered zero-crossing). The velocity zero-crossing frequency is statistically analyzed within each time window: The number of zero-crossings within this window is counted and divided by the time window coverage duration to obtain the frequency of that window, thus obtaining the velocity zero-crossing frequency sequence.

[0056] Frequency matching is calculated window-by-window, aligned to the time window. For each time window, the difference between the curvature zero-crossing frequency and the velocity zero-crossing frequency is first calculated, and the difference is taken as the absolute difference between the two. Then, the difference is normalized and mapped to a matching score: the difference is divided by the acceptable difference upper limit to obtain the normalized difference. When the normalized difference is 0, the matching score is 1; when the normalized difference reaches 1, the matching score is 0; and when the normalized difference exceeds 1, the matching score is 0 to avoid negative values. The acceptable difference upper limit is determined during the modeling phase: the two types of frequency differences are calculated window-by-window for strongly coupled trajectory segments, and the maximum difference is taken and multiplied by 1.1 as the upper limit (the multiplier of 1.1 is used to cover operating condition fluctuations and ensure that strongly coupled samples are not misjudged as mismatches). Strongly coupled trajectory segments are screened by the previous coupling judgment results and included in the statistics after manual sampling to confirm label consistency.

[0057] The distribution trends of the bending segment length and acceleration amplitude are calculated within the same time window, and the 64-point time window and 16-point step size are shared with the frequency matching degree to ensure alignment. Bending segment division is based on the curvature amplitude threshold: for each time window, the bending threshold is set to the lower limit of curvature noise multiplied by 3 (the multiplier of 3 is determined by the interval distinguishing between straight-line noise and actual steering: the distribution of the lower limit of curvature noise in straight-line segments and the curvature amplitude of steering segments are statistically analyzed on the modeling data; the multiplier of 3 minimizes the overlap between the two distributions and maintains a steering segment detection rate of no less than 0.9). The interval of sampling points within the window that continuously satisfy the curvature amplitude not lower than the bending threshold is marked as a bending segment. The bending segment length is obtained by accumulating the spatial distances of adjacent sampling points within the interval: for each pair of adjacent sampling points within the bending segment, a three-dimensional coordinate difference vector is calculated. The square root of the sum of the squares of each component of the difference vector is used to obtain the distance of that segment. Then, the distances of all adjacent points covered by the bending segment are accumulated segment by segment to obtain the bending segment length. A set of bending segment lengths is formed within the time window. The acceleration amplitude set is directly composed of the acceleration amplitudes of the 64 sampling points within the window. Subsequently, the distribution of bending segment length and acceleration amplitude are constructed separately. The distribution uses eight equally wide intervals for statistical analysis (the number of intervals, 8, is determined by both resolution and computational cost: after evaluating the stability and computational delay of the distribution similarity discrimination using candidate interval numbers of 6, 8, and 10, the number of intervals, when set to 8, achieves optimal sensitivity and stability to changes in operating conditions while minimizing computational delay). The boundaries of the bending segment length intervals are determined by the strongly coupled data from the modeling phase: the minimum and maximum values ​​of the bending segment lengths are statistically analyzed, and the range from the minimum to the maximum value is divided into eight equal segments to form the boundary; the boundaries of the acceleration amplitude intervals are formed using the same method, dividing the minimum and maximum acceleration amplitude values ​​into eight equal segments. For the bending segment length set, each bending segment length falls into its corresponding interval, and the interval is counted, resulting in an 8-dimensional counting vector; for the acceleration amplitude set, the acceleration amplitude of each sampling point falls into its corresponding interval, and the interval is counted, resulting in an 8-dimensional counting vector. The two sets of count vectors are normalized into probability vectors. Normalization is achieved by dividing the count for each interval by the sum of the counts in that vector, yielding the proportion of each interval and ensuring comparability between the two distributions. The distribution matching score is obtained through probability difference mapping: first, the absolute difference between the two probabilities is calculated for each of the eight intervals, and these differences are accumulated to obtain the difference quantity; the theoretical upper limit of the difference quantity is 2, so the difference quantity is divided by 2 to obtain the normalized difference; then, 1 is subtracted from the normalized difference to obtain the matching score, ensuring a score of 1 for perfect agreement and 0 for perfect disagreement. The distribution matching trend is constituted by adjacent changes in the matching score sequence: for each time window, the current matching score is recorded, and the change is calculated by subtracting the matching score of the previous time window from the current matching score. A positive change indicates increased matching, and a negative change indicates decreased matching; the change for the first time window is set to 0 to maintain a consistent sequence length.

[0058] The temporal feature sequence is obtained by concatenating the frequency matching degree and distribution matching information of each time window according to fixed fields. One temporal feature vector is output for each time window, containing four components: current frequency matching degree, frequency matching degree change, current distribution matching score, and distribution matching score change. To eliminate the impact of amplitude differences in different trajectory segments on network state updates, standardization is performed on the four components. The standardization parameters are taken from the statistics of the modeling phase: the mean and standard deviation are calculated for each component. The mean is obtained by summing the training samples item by item and dividing by the number of samples. The standard deviation is obtained by summing the squares of each item after removing the mean, dividing by the number of samples minus 1, and then taking the square root. The lower limit of the standard deviation is also set to 0.001 (the determination of the value 0.001 is based on the same criteria as before, used to suppress standardization amplification noise). Standardization is performed item by item for each component, subtracting the corresponding mean from each component and then dividing by the corresponding standard deviation to obtain the standardized temporal feature sequence.

[0059] The standardized temporal feature sequence enters the Long Short-Term Memory (LSTM) network. The network processes each feature vector sequentially according to time windows and updates the memory state and hidden state. For any time window, the network first concatenates the current 4-dimensional feature vector with the hidden state vector of the previous time window to form a joint input. The length of the hidden state vector is set to 128 (the value 128 is jointly determined by the dynamic pattern expressiveness and inference latency: after evaluating the trajectory pattern discrimination accuracy and single-segment inference latency, the candidate lengths of 64, 96, and 128 show that when the length is 128, the discrimination accuracy improvement reaches a stable range and the inference latency does not exceed the budget). Based on the joint input, the gate values ​​and candidate values ​​of the forget gate, input gate, candidate memory, and output gate are calculated respectively: each component of each type of gate value is obtained by multiplying the joint input component with the corresponding weight item by item and then summing them, and then superimposing the bias to obtain the pre-gate activation value; the pre-gate activation value is transformed by a nonlinear function with a limited output range of 0 to 1 to obtain the gate value, which is used to control the information retention and introduction intensity. Each component of a candidate memory is obtained by multiplying and adding the joint input and corresponding weights, and then superimposing the bias to obtain the candidate pre-activation value. The candidate pre-activation value is then transformed by a nonlinear function with an output range of -1 to +1 to obtain the candidate memory value, which is used to provide new memory content. Memory state updates are performed component by component: each component of the memory state of the previous time window is first multiplied by the corresponding forget gate value to obtain the retained part, then the candidate memory value is multiplied by the input gate value component by component to obtain the added part, and finally the retained part and the added part are added component by component to obtain the current memory state. Hidden state updates are performed component by component: the current memory state is first transformed by a nonlinear function with an output range of -1 to +1 to obtain the memory output, and then multiplied by the output gate value component by component to obtain the current hidden state. The aforementioned weights and biases are determined through supervised learning during the modeling phase. The supervision objective is defined as making the time-series dependence of the hidden state on the frequency matching degree and distribution trend form a stable representation and improving the discriminativeness of the trajectory dynamics pattern. The training data consists of labeled multi-condition trajectory segments. The optimization process selects the final parameters based on minimizing the validation set error. The learning rate is set to 0.001 (the value of 0.001 is determined by the combined constraints of the loss descent rate and oscillation amplitude. When 0.001 is selected from the candidates 0.0005, 0.001, and 0.002, convergence is stable and the validation error is the lowest). The weight decay coefficient is set to 0.0001 (the value of 0.0001 is determined by the overfitting suppression requirement, so that the difference between the training error and the validation error does not exceed 0.02 and the pattern discriminativeness does not decrease).

[0060] The trajectory dynamics pattern vector is obtained by mapping the hidden state vector output by the Long Short-Term Memory network. The mapping takes the hidden state vector corresponding to the last time window as the dynamic representation of the entire segment. The mapping process adopts a fixed linear transformation: for the first component of the pattern vector, the 128 components of the hidden state vector are multiplied by their corresponding weights one by one and then accumulated, and the bias is added to obtain the output; the multiplication and addition process is repeated for the second to the 128th components to obtain a trajectory dynamics pattern vector of length 128 (the length 128 is determined by the pattern capacity requirements of the subsequent coupling verification and trajectory optimization module: after evaluating the pattern clustering separation and subsequent optimization convergence stability, the candidate lengths of 64, 96, and 128 are selected. When the length is 128, the separation improvement reaches the stable range and the mapping calculation delay meets the budget). To ensure comparability of the pattern vector across different trajectory segments, amplitude shaping is performed on the pattern vector. Amplitude shaping employs component-by-component pruning and normalization: for each component, if it exceeds the upper bound statistically determined during the modeling phase, the upper bound is used; if it falls below the lower bound, the lower bound is used. The upper and lower bounds are obtained by multiplying the maximum and minimum values ​​of the training sample components by 1.1 (the multiplier of 1.1 is used to cover fluctuations in operating conditions and suppress occasional outliers). Subsequently, standardization is performed based on the mean and standard deviation statistically determined during the modeling phase, with the lower limit of the standard deviation still set at 0.001 (the determination of the value 0.001 is based on the same criteria as before). The trajectory dynamics pattern vector obtained through the above processing is then bound to the trajectory segment corresponding to the time index.

[0061] S5 includes acquiring the trajectory dynamics mode vector and geometric shape feature vector, calculating the path fractal dimension and dynamic energy dissipation rate respectively; calculating the mutual information value of the path fractal dimension and dynamic energy dissipation rate as the correspondence score; if the correspondence score is lower than a preset threshold, updating the control point position along the gradient vector that maximizes the correspondence score to reconstruct the corrected industrial robot trajectory sample point sequence; performing third-order differentiation on the corrected industrial robot trajectory sample point sequence to obtain the jerk sequence; if the jerk sequence satisfies the smoothness limit value, it is determined that the adjusted sequence satisfies the predefined kinematic continuity constraint.

[0062] In this embodiment, the trajectory dynamics mode vector and geometric feature vector are obtained simultaneously for the same trajectory segment, and the sequence of trajectory sample points for that trajectory segment is used as the calculation benchmark. The trajectory segment length is set to 256 sampling points (256 is determined by the controller interpolation cycle of 0.004 seconds and the calculation delay budget of a single planning cycle; using 256 points covers the continuous interval of steering and acceleration / deceleration combined actions, and the single-segment processing delay is consistently lower than the planning cycle budget). The sampling interval is set to 0.004 seconds (0.004 is determined by reading the controller interpolation cycle and serves as a unified time benchmark for velocity, acceleration, and higher-order differences). The trajectory sample point sequence contains the three-dimensional position, three-dimensional velocity, and three-dimensional acceleration of each sampling point. The three-dimensional velocity is obtained by dividing the position difference of two adjacent sampling points by 0.004 seconds, and the three-dimensional acceleration is obtained by dividing the velocity difference of two adjacent sampling points by 0.004 seconds. Boundary sampling points use a weighted combination of the differences of two adjacent segments to avoid amplifying noise from the boundary difference. The weights are 0.6 and 0.4 (0.6 and 0.4 are jointly determined by boundary difference noise suppression and velocity continuity constraints to ensure that the velocity and acceleration of the first point do not jump abruptly and do not introduce leading information).

[0063] The path fractal dimension is calculated based on multi-scale coverage statistics of the trajectory segment position sequence. First, the minimum bounding box of the trajectory segment in the workspace is calculated. The boundaries of the bounding box in the three coordinate directions are obtained by comparing the minimum and maximum values ​​of all position sample points in the segment in the corresponding directions. Then, eight grid scales are set (eight is determined by the stability of dimension estimation and the budget for single-segment computation delay; fewer than eight scales will lead to increased dimension fluctuations, and more than eight scales will lead to grid indexing and occupancy statistics delays exceeding the budget). The maximum scale is half the length of the longest side of the bounding box (half is taken to avoid the number of occupied grids degenerating to a constant due to excessively large scales), and the minimum scale is the length of the longest side of the bounding box divided by 64 (64 is determined by the sampling space resolution and detail preservation requirements to ensure that the minimum grid side length is not less than the average distance between adjacent sampling points, thereby avoiding single-point noise dominating the occupancy count). For each scale, the bounding box is divided into three-dimensional grid cells according to that scale, and grid localization is performed on sample points at each location of the trajectory segment: the coordinates of the sample points are subtracted from the minimum boundary of the bounding box in each direction, then divided by the grid edge length in each direction and rounded down to obtain the grid index; after indexing all sample points, the number of grid cells containing at least one sample point is counted as the coverage quantity for that scale. After completing the coverage quantity for 8 scales, in order to reduce the influence of endpoint saturation and endpoint noise, one set of data each from the largest and smallest scales is removed, and only the middle 6 sets of data are retained for dimension estimation (the removal rule is determined by the statistical characteristic that the coverage quantity at the endpoint scale is more susceptible to saturation and single-point noise, and the variance of the dimension result on the historical trajectory segment is significantly reduced after removal). Dimension estimation is achieved using a linear fitting slope: First, six sets of scales and coverage quantities are mapped to a logarithmic scale. The mapping process obtains corresponding values ​​through a pre-set lookup table, avoiding complex numerical calculations on-site. Then, least squares fitting is performed on the six sets of logarithmic scales and logarithmic coverage quantities. The fitting calculation is completed by point-by-point accumulation, including the logarithmic scale mean, the logarithmic coverage quantity mean, the cumulative sum of deviations, and the cumulative sum of squared deviations. Finally, the slope is obtained by dividing the cumulative sum of deviations by the cumulative sum of squared deviations. This slope is used as the path fractal dimension output, so that the dimension value is higher when the path shape is more complex.

[0064] The dynamic energy dissipation rate is a segment-level index based on the velocity and acceleration amplitudes of the trajectory segment, and is modulated by the trajectory dynamics mode vector. The velocity amplitude is obtained by summing the squares of the three velocity components at each sampling point and then taking the square root; the acceleration amplitude is obtained by summing the squares of the three acceleration components at each sampling point and then taking the square root. For each sampling point, the velocity amplitude and acceleration amplitude are multiplied to obtain a power intensity sample. Then, the power intensity samples for all points in the trajectory segment are summed point by point, and the result is divided by the number of sampling points (256) to obtain the segment-level power intensity. The mean value is taken to suppress the interference of instantaneous spikes on the segment-level dissipation judgment, ensuring that the index reflects the overall dissipation level. The trajectory dynamics mode vector is mapped to obtain modulation coefficients: each component of the mode vector is multiplied by its corresponding weight, the results are accumulated, and a bias is added to obtain the original coefficients. These original coefficients are then cropped to a range of 0.5 to 2.0 (0.5 and 2.0 are determined by the energy dissipation rate; values ​​below 0.5 will compress the differences in high-energy-consuming segments, leading to decreased discrimination, while values ​​above 2.0 will amplify noise, causing increased score fluctuations; this range is selected in the correlation evaluation between historical trajectory segments and the mean square value of the driving current, maximizing correlation and minimizing fluctuations). The segment-level power intensity is multiplied by the modulation coefficients to obtain the dynamic energy dissipation rate output.

[0065] The correlation score between path fractal dimension and dynamic energy dissipation rate is represented by mutual information, which is calculated using sliding trajectory segment sample pairs as input. The same trajectory is divided into a continuous sequence of trajectory segments by a sliding window. The window length is 256 sampling points, and the step size is 16 sampling points (16 is determined by the number of statistical samples and the smoothness of adjacent windows; with a step size of 16, adjacent windows overlap by 240 sampling points, resulting in a smooth mutual information sequence and forming enough sample pairs within a fixed trajectory length for stable estimation). For each window, a path fractal dimension and a dynamic energy dissipation rate are calculated, forming a set of sample pairs. Subsequently, binning discretization is performed on both indices, with a bin size of 16 (16 is determined by the mutual information estimation accuracy and the risk of joint frequency sparsity; too small a bin size will lose dependency structure, while too large a bin size will lead to joint count sparsity and introduce estimation bias; 16 minimizes mutual information fluctuations within the historical sample size). The binning boundaries are obtained by dividing the minimum and maximum value intervals of strongly coupled trajectory segments with the same width (using strongly coupled trajectory segments to cover typical operating conditions, and dividing with the same width to ensure the uniqueness and reproducibility of binning rules at the deployment end). After discretization, the frequency of occurrence of joint bin pairs is counted to form a joint frequency table, and the bin frequencies of each of the two types of indicators are counted separately to form an edge frequency table. The joint frequency and edge frequency are divided by the total number of sample pairs to obtain the joint probability and edge probability, respectively. Mutual information accumulation is performed item by item for each bin pair: for each bin pair with a joint probability greater than 0, the ratio of the product of the joint probability and the corresponding edge probability of the bin pair is calculated, and then the natural logarithm of the ratio is taken. Subsequently, the natural logarithm is multiplied by the joint probability of the bin pair and accumulated to the total mutual information sum; bin pairs with a joint probability of 0 are not included in the accumulation to avoid introducing numerical anomalies due to invalid terms. The sum of mutual information is linearly normalized to obtain the corresponding score: the upper bound of normalization is the maximum mutual information of the strongly coupled trajectory segment multiplied by 1.1 (1.1 is determined by the upper bound of the operating condition fluctuation, which is used to avoid the strongly coupled trajectory segment being suppressed by the field noise and causing false triggering correction), and the lower bound of normalization is 0. After normalization, the score is clipped to the range of 0 to 1, which is used as the corresponding score of the path fractal dimension and the dynamic energy dissipation rate.

[0066] Control point updates are triggered when the corresponding relationship score falls below a preset threshold of 0.6 (0.6 is determined by a combination of safety and correctability constraints: the historical trajectory segment is scanned in increments of 0.05 steps from 0.4 to 0.8; the proportion of tracks that, after triggering, have no more than 30 iterations and whose scores improve to 0.75 or higher is no less than 0.9; and the geometric offset of the original qualified trajectory does not exceed 2 mm; the final threshold of 0.6 satisfies the above constraints). Control points are obtained using cubic spline parameterization, and the number of control points is 32 (32 is determined by the local plasticity of the trajectory and the scale of the optimization variables; if the number is less than 32, the bending change required for local obstacle avoidance correction is insufficient; if the number is greater than 32, the gradient estimation and reconstruction calculation delay exceeds the budget). The control point update direction is the gradient direction that maximizes the correspondence score. The gradient is obtained component-by-component using a finite difference approximation: For each coordinate component of each control point, first increase by 0.1 mm along the positive direction of the coordinate and reconstruct the trajectory sample point sequence accordingly, then calculate the perturbation-adjusted correspondence score; then decrease by 0.1 mm along the negative direction of the coordinate and reconstruct the trajectory sample point sequence, then calculate the perturbation-adjusted correspondence score again; divide the difference between the two scores by 0.2 mm to obtain the gradient component of that coordinate component (0.1 mm is determined by both control accuracy and differential stability; when the perturbation is less than 0.1 mm, the score change is drowned out by numerical noise, and when the perturbation is greater than 0.1 mm, the gradient approximates the local linear range). The gradient components of the three coordinate components of all control points are concatenated in the order of the control points to form the gradient direction. The control point position is updated according to the gradient direction, with an update step size of 0.5 mm (0.5 mm is determined by both convergence speed and trajectory morphology stability; a step size less than 0.5 mm leads to more than 30 iterations, and a step size greater than 0.5 mm causes local morphological abrupt changes and increases the risk of jerk exceeding limits). After each update, spatial constraint correction is immediately performed: the control point position must not enter within a 20mm safety distance of obstacles (20mm is determined by the sum of the maximum radius of the end effector, the upper limit of the repeatability error, and the upper limit of the sensor positioning error, plus 5mm; the 5mm is used to cover instantaneous offsets caused by load changes). If it enters this range, the control point is pushed back to the boundary position 20mm away along the outer normal direction of the obstacle to ensure that the control point update does not introduce collision risks. The updated control points are used to reconstruct the corrected trajectory sample point sequence. The reconstruction uses cubic spline interpolation: for adjacent control point intervals, the position at each sampling moment within the interval is calculated based on the continuity constraint of the interval endpoint position and the first derivative of the endpoint, ensuring that the position, velocity, and acceleration are continuous at the interval splicing, thus forming the corrected industrial robot trajectory sample point sequence.

[0067] The corrected sequence of industrial robot trajectory sample points must satisfy predefined kinematic continuity constraints, with continuity determined by the smoothness limit of the jerk sequence. First, the velocity sequence is calculated from the corrected position sequence with a 0.004-second difference. Then, the acceleration sequence is calculated from the velocity sequence with a 0.004-second difference. Finally, the jerk sequence is obtained from the acceleration sequence with a 0.004-second difference. The difference process is completed by subtracting each component from adjacent sampling points and dividing by 0.004 seconds. Boundary points are obtained by weighted combinations of two adjacent difference segments, using weights of 0.6 and 0.4 to suppress boundary noise. The jerk amplitude is calculated for each sampling point, obtained by summing the squares of the three jerk components and taking the square root. The smoothness limit is the maximum allowable jerk amplitude of the controller multiplied by 0.9 (0.9 is determined by a safety margin to cover instantaneous overshoot caused by load and friction changes, avoiding triggering controller protection during actual execution). When the acceleration amplitude of all sampling points within the trajectory segment does not exceed the limit value, the corrected sequence is deemed to satisfy the kinematic continuity constraint; if any sampling point exceeds the limit, it is deemed not to satisfy the constraint and the next control point update iteration begins. The upper limit of the number of iterations is 30 (30 is jointly determined by the total time delay budget of the planning cycle and the computational cost of a single gradient estimation, reconstruction, and mutual information calculation; when the upper limit is 30, the total processing delay is still lower than the planning cycle budget in the worst case). If the smoothness limit value is not satisfied even after reaching the upper limit, the set of control points with the highest score in the last time the corresponding relationship is satisfied is maintained, and the control point update step size is reduced from 0.5 mm to 0.25 mm before continuing execution until the smoothness limit value is satisfied (0.25 mm is determined by the half-step strategy, used to suppress the sudden increase of higher-order changes without changing the gradient direction).

[0068] S6 includes acquiring the trajectory sequence of the industrial robot and the obstacle model that satisfy the kinematic continuity constraints, identifying potential obstacle regions in the trajectory sequence located near the obstacle model; calculating the synchronization degree based on the local curvature changes and acceleration peaks in the potential obstacle regions to determine the interference probability; if the interference probability is higher than a preset threshold, using a dynamic programming algorithm to search for the optimal state transition path passing through candidate nodes outside the potential obstacle regions to obtain the preferred obstacle avoidance path.

[0069] In this embodiment, the input consists of an industrial robot trajectory sequence and an obstacle model that satisfy kinematic continuity constraints. The trajectory sequence is arranged chronologically, and each sampling point contains three-dimensional position, three-dimensional velocity, and three-dimensional acceleration information. The sampling interval is 0.004 seconds (0.004 seconds is determined by the controller's interpolation cycle and serves as a unified time reference for distance calculation, curvature calculation, peak detection, event alignment, and path resampling). The obstacle model is represented by a three-dimensional surface mesh or voxel set. The obstacle model coordinate system and the trajectory sequence coordinate system are the robot's base coordinate system (coordinate system one is confirmed by the controller's currently effective coordinate configuration, and the confirmation method is that the coordinate transformation hysteresis of the same spatial point under the two coordinate configurations is no greater than 0.05 mm, and 0.05 mm is determined by superimposing the coordinate transformation numerical accuracy and the repeatability of on-site calibration).

[0070] First, potential obstacle regions are identified. A safety expansion is performed on the obstacle model to obtain a safety envelope. The safety expansion distance is 20 mm (this 20 mm is determined by adding the maximum radius of the end-effector, the upper bound of the repeatability error, and the upper bound of the environmental perception positioning error, plus 5 mm; the 5 mm is used to cover instantaneous deviations caused by load changes, ensuring a deterministic safety gap when the trajectory approaches the obstacle). When the obstacle model is represented as a surface mesh, safety expansion generates offset patches by shifting the outward normal direction of each triangular facet of the mesh by 20 mm. Then, gap repair and self-crossing correction are performed on the set of offset patches to form a closed envelope surface. Gap repair is completed by pairing the nearest points of adjacent patch boundaries and inserting connecting patches. Self-crossing correction is completed by deleting internal patches after distinguishing between the inside and outside of the envelope surface. When the obstacle model is expressed as a set of voxels, safety expansion is achieved by expanding the three-dimensional neighborhood. The number of expansion layers is the safety expansion distance divided by the voxel side length and then rounded up. Rounding up is used to ensure that the expansion distance is not less than 20 mm. The expansion is implemented by marking the newly added occupying voxel in its three-dimensional neighborhood for each occupying voxel until the expansion layer is completed.

[0071] After the safety envelope is established, the shortest distance from each trajectory point to the safety envelope is calculated for each point in the trajectory sequence. When using surface mesh representation, a spatial partition index is first established for the safety envelope. The side length of the partition mesh is 50 mm (50 mm is determined by the upper bound of the obstacle density and the real-time calculation budget; when the side length is 50 mm, the number of candidate patches for each query is stable and controlled, and the nearest patch is not missed). The partition index is established by mapping the bounding box of the envelope patch to the partition mesh index and registering the patch number in the corresponding cell. For each trajectory point, first locate its corresponding grid cell, then take the candidate faces registered in that grid cell and its adjacent grid cells, calculate the shortest distance from the point to the triangular face for each face, and take the minimum value as the shortest distance of the trajectory point; the calculation of the shortest distance from the point to the triangular face is completed by projection discrimination. First, project the point onto the plane where the face is located. If the projected point falls inside the triangle, the shortest distance is the perpendicular distance from the point to the plane, which is calculated from the normal projection length from the point to the plane; if the projected point falls outside the triangle, the shortest distance is the minimum value of the shortest distance from the point to the three side line segments of the triangle. The shortest distance of the line segments is obtained by calculating the distance from the point to the endpoint of the line segment or the point to the projection point of the line segment after truncation by the line segment parameters. When using voxel sets for representation, the coordinates of each trajectory point are converted into a voxel index, and the occupancy of neighboring voxels centered on that index is queried. The neighborhood radius is converted into the number of voxel grids by the safe expansion distance and then rounded up. If there is an occupied voxel in the neighborhood, the nearest distance is approximated as the distance from the trajectory point to the center of the nearest occupied voxel. The upper bound of the error of this approximation is determined by half the length of the voxel diagonal, and a margin of 20 mm has been reserved for coverage.

[0072] Track points with a nearest distance of no more than 40 mm are marked as near-obstacle points (40 mm is determined by a safety expansion distance of 20 mm plus an early warning distance of 20 mm; the early warning distance is used to cover the spatial offset caused by subsequent detour path generation and resampling). Near-obstacle points are merged according to their time indices, with the merging rule being that the time index interval between adjacent obstacle points does not exceed 3 sampling points (these 3 sampling points are jointly determined by sampling jitter and the upper bound of short-term fluctuations in distance calculation, to avoid the near-obstacle region being fragmented into multiple fragmented intervals by fluctuations). After merging, several potential obstacle regions are formed, each region containing a starting sampling point index, an ending sampling point index, and the corresponding trajectory segment position, velocity, and acceleration sequence.

[0073] Subsequently, the synchronization degree between local curvature changes and peak acceleration is calculated within each potential obstacle region, and the interference probability is determined based on the synchronization degree. The local curvature change sequence is generated point-by-point based on the position sequence within the region. For the second to the second-to-last sampling point within the region, a three-point neighborhood is formed by taking the previous point, the current point, and the next point. First, the displacement vectors of the preceding and following segments are calculated. Each component of the displacement vector is obtained by subtracting the corresponding coordinates, and the displacement length is obtained by taking the square root of the sum of the squares of each component. Then, each component of the displacement vector is divided by the displacement length to obtain the unit direction vector. Subsequently, the turning amplitude of the unit direction vectors of the preceding and following segments is calculated. The turning amplitude is obtained by converting the dot product of the two vectors, and the dot product is obtained by multiplying the corresponding components and summing them one by one. The turning amplitude is normalized according to the neighborhood average step size to obtain the bending strength. The neighborhood average step size is the arithmetic mean of the preceding and following displacement lengths. The arithmetic mean is used to suppress the amplification of bending strength caused by unilateral step size anomalies. The local curvature variation amplitude is obtained by the adjacent differences of the bending intensity. For each sampling point, the variation amplitude is obtained by subtracting the bending intensity of the previous point from the current bending intensity and arranged in time to form a sequence. To suppress the influence of single-point impulse noise on event detection, median filtering is performed on the variation amplitude sequence. The filtering window length is 3 sampling points (3 is determined by the statistical characteristic that impulse noise often manifests as a single-point spike, while constraining the event position offset to not exceed 1 sampling point).

[0074] The acceleration peak sequence is generated point-by-point based on the acceleration amplitude sequence within the region. The acceleration amplitude is calculated for each sampling point, obtained by summing the square roots of the squares of the three-dimensional acceleration components. Peak detection uses a sliding window with a window length of 5 sampling points (5 is determined by the upper bound of the maximum event offset after converting servo response hysteresis and sampling interval, ensuring that the true peak falls within the window coverage). For each window, if the amplitude at the center point is strictly greater than the amplitudes of the other 4 points within the window, the center point is marked as the peak point and its amplitude is output; the remaining points are output as 0, forming the peak sequence. To avoid excessive peak density due to noise, a minimum peak interval constraint is introduced, with a minimum peak interval of 3 sampling points (3 is determined by the shortest significant acceleration change interval achievable by the mechanical system; intervals less than 3 will result in physically meaningless repetitive peaks).

[0075] Synchronization is calculated using an event matching method. First, curvature events are detected in the local curvature change amplitude sequence. The detection window length is set to 5 sampling points (5 is consistent with the peak detection window to maintain consistent temporal resolution for both types of events). Events with a center point change amplitude significantly greater than the other four points within the window are marked as positive curvature events, while those with a center point change amplitude significantly less than the other four points are marked as negative curvature events. Then, using each curvature event as a baseline, acceleration peak points are searched within a range of two sampling points before and after it (these two sampling points are determined by the maximum alignment offset from historical data statistics; the maximum offset is set to two sampling points, and the range is consistent with the maximum offset to cover hysteresis and jitter). If a peak point exists within the search range, the peak point with the smallest time offset is selected as the matching point, and a matching score is given according to the time offset: 1.0 is awarded when the offset is 0 sampling points (1.0 is the full score, used to indicate strict alignment); 0.7 is awarded when the offset is 1 sampling point (0.7 is determined by the alignment error distribution of strongly coupled samples, with a high proportion of samples having an offset of 1 and needing to maintain a significant contribution); 0.4 is awarded when the offset is 2 sampling points (0.4 is determined by the decreasing proportion of samples having an offset of 2 and needing to increase the distinguishability). If no peak point exists within the search range, the score for that curvature event is 0. The synchronization degree value is obtained by summing the scores of all curvature events within the region and dividing by the number of curvature events; when the number of curvature events is 0, the synchronization degree value is 0 to avoid a denominator of 0 and to avoid artificially high synchronization in areas without events.

[0076] The interference probability is obtained by mapping the synchronization level value and the obstacle proximity level together. The obstacle proximity level consists of the minimum and average values ​​of the nearest distance within the area: the minimum value is obtained by comparing points one by one, and the average value is obtained by summing points one by one and dividing by the number of sampling points in the area. The interference probability is obtained by inputting the synchronization level value, the minimum nearest distance value, and the average nearest distance value into the interference probability mapping table. The interference probability mapping table uses discrete tiered storage, with 5 tiers for synchronization level, 5 tiers for minimum nearest distance value, and 5 tiers for average nearest distance value, for a total of 125 items; the tier boundaries are jointly determined by historical collision alarm data and no alarm data, and the determination method is to ensure that the alarm detection rate is not less than 0.9 under the constraint that the false alarm rate is not higher than 0.1, and that the interference probability increases monotonically with the increase of synchronization level and increases monotonically with the decrease of nearest distance (the upper limit of false alarm rate of 0.1 and the lower limit of detection rate of 0.9 are determined by obstacle avoidance safety requirements). During runtime, the table entry is directly located and the interference probability is read based on the index of the tiered interval into which the input value falls. The obtained interference probability is then clipped to the range of 0 to 1, with the clipping rule being 0 for values ​​less than 0 and 1 for values ​​greater than 1.

[0077] When the interference probability exceeds a preset threshold of 0.6, path replanning is triggered (0.6 is determined by a combination of safety and path disturbance magnitude constraints: on the validation data, the threshold is iterated from 0.4 to 0.8 in increments of 0.05, and thresholds that meet the following criteria are selected: obstacle avoidance success rate not less than 0.98, path length increment not exceeding 0.15 of the original path length, and planning calculation delay not exceeding the planning cycle budget). Path replanning uses dynamic programming to search for the optimal state transition path outside the potential obstacle area. The outer boundary is obtained by expanding the safety envelope by 30 mm (30 mm is determined by the superposition of the detour channel reserve and the maximum spatial drift after trajectory correction, ensuring that the search space includes feasible detour corridors and controlling the search scale). A candidate node set is constructed within the outer boundary. Node generation adopts hierarchical sampling: within the time index range of the potential obstacle area along the original trajectory, one reference point is selected every 8 sampling points (8 is determined by the joint constraints of node density and the number of dynamic programming states; when 8 is selected, the number of states is controlled and the turning granularity of the bypass path satisfies the curvature constraint); with each reference point as the center, the local tangential direction of that point is calculated. The tangential direction is obtained by normalizing the position difference between that point and the previous point. The normalization process is to divide each component of the difference vector by the difference length, and the difference length is obtained by summing the squares of each component and taking the square root. Then, a normal plane orthogonal to the tangential direction is constructed. Within this plane, three lateral nodes are generated along the main offset direction away from the obstacle. The lateral offset distances are 10 mm, 20 mm, and 30 mm (10, 20, and 30 are determined by the safety envelope expansion gradient and trajectory smoothness, with 10 corresponding to micro-detour and 30 corresponding to strong detour, covering typical obstacle avoidance offset amplitudes). At the same time, nodes at the same distance are generated on both sides of the main offset direction, each deflected by 30 degrees. The 30 degrees are determined by the detour channel width and the upper limit of curvature, ensuring that the nodes cover the outer channel of the obstacle without introducing sharp backtracking. Feasibility screening is performed on the generated nodes: the shortest distance from the node to the safety envelope must be no less than 20 mm (20 mm is taken from the safety expansion distance to ensure that the node falls outside the defined safety gap); the local steering intensity formed by the connection between adjacent reference layer nodes at the node must not exceed the upper limit of curvature. The upper limit of curvature is determined by converting the minimum turning radius allowed by the controller with the upper limit of the planned speed in this area. The conversion is based on the fact that the smaller the turning radius under the upper limit of the planned speed, the greater the steering intensity. The maximum steering intensity corresponding to the minimum turning radius allowed by the controller is taken as the upper limit; nodes that do not meet the screening conditions are directly removed to avoid entering the dynamic programming state set.

[0078] In dynamic programming, a state is defined as "a candidate node in a certain reference layer," and a state transition is defined as a connection segment formed by connecting a node in the previous reference layer to a node in the next reference layer. Constraints are evaluated and transition costs are calculated for each candidate connection segment. Collision constraint determination uses connection segment sampling detection, with 5 sampling points (5 is determined by both the cost integral approximation accuracy and computational complexity; 5 points ensure stable detection of crossing risks in narrow passages while keeping computational complexity under control). Five sampling points are taken at equal intervals on the connection segment, with the interval value obtained by dividing the connection segment length by 6 to ensure that sampling points cover the endpoints and the middle. The nearest distance to the safety envelope is calculated for each sampling point, using the same method as for near-obstacle recognition. If the nearest distance to any sampling point is less than 20 mm, the connection segment is deemed unusable. The continuity constraint is determined based on the upper limit of tangential change. Tangential change at the start of the connecting segment is calculated by converting the angle between the direction of the previous segment and the direction of the connecting segment, and at the end of the connecting segment, it is calculated by converting the angle between the direction of the connecting segment and the direction of the next segment. The dot product required for angle conversion is obtained by multiplying the corresponding components of the direction vector and then summing them one by one. The length of the direction vector is obtained by taking the square root of the sum of the squares of the corresponding components. The upper limit of tangential change is derived from the preceding kinematic continuity constraints and the jerk limit, and verified on the validation trajectory to ensure that no high-order abrupt changes occur at the connection. The risk constraint is determined based on the cumulative value of the interference probability along the connecting segment. The interference probability along the segment is obtained by querying the interference probability mapping table for each of the five sampling points. The cumulative value is the arithmetic mean of the five probabilities multiplied by the length of the connecting segment. The length of the connecting segment is obtained by taking the square root of the sum of the squares of the differences between the endpoint coordinates. The upper limit of risk is taken as the lower bound of the risk level corresponding to the trigger threshold of 0.6, ensuring that the overall risk of the bypass segment is lower than the risk of the trigger segment.

[0079] The transfer cost consists of three parts: path length cost, risk cost, and smoothing cost, which are summed in weights. The path length cost is the length of the connecting segment. The risk cost is the cumulative value of the aforementioned risks. The smoothing cost is the sum of the tangential changes at the start and end points of the connecting segment. The weights for the three components are 0.5, 0.3, and 0.2 (these values ​​are determined by an optimal trade-off under multi-objective constraints: the obstacle avoidance success rate is the primary constraint; the path length increment is minimized while ensuring the success rate is not less than 0.98, and the smoothing cost is limited to the controller's executable range; the weight combinations are iterated through on the validation data, and the combination that satisfies the constraints and minimizes the total cost is selected). Dynamic programming employs a layer-by-layer cumulative minimum cost approach: In the first reference layer, the initial cost of entering a node from the original trajectory is calculated for each node and recorded. In the next reference layer, for each node, all nodes in the previous reference layer are traversed, and the sum of the minimum cumulative cost of the previous node and the transition cost between the two nodes is calculated. The minimum cost is taken as the minimum cumulative cost of that node, and the index of the predecessor node that generated the minimum cost is recorded. This process is repeated in the order of the reference layers until the final reference layer. After completing the calculation in the final reference layer, the node corresponding to the minimum cumulative cost is selected as the endpoint node. The node sequence is obtained by backtracking along the predecessor node index, and the node sequence is connected to form the obstacle avoidance optimal path.

[0080] Before outputting the optimal obstacle avoidance path, smoothing and resampling processes are performed to maintain continuity with the original trajectory. Smoothing employs cubic spline fitting, with the input being a backtracked sequence of nodes. Fitting constraints include the starting position aligning with the original trajectory's entry point, the ending position aligning with the original trajectory's exit point, the starting tangent aligning with the original trajectory's entry tangent, and the ending tangent aligning with the original trajectory's exit tangent. The tangent is obtained by normalizing the differences between adjacent sampling points, with each component of the difference vector divided by the difference length. After fitting, resampling is performed at 0.004-second intervals. Resampling generates time-step position points by advancing the spline curve by arc length, and generates velocity and acceleration based on the position point differences, ensuring the new path sequence aligns with the controller's interpolation timing, thus forming the final optimal obstacle avoidance path sequence.

[0081] S7 includes calculating the spatial overlap density between curvature extrema and velocity mutation points in the obstacle avoidance optimization path, and extracting local trajectory segments based on the spatial overlap density; obtaining the geometric self-similarity and acceleration sequence autocorrelation peak of the local trajectory segments, and constructing a multidimensional feature tensor based on the geometric self-similarity and acceleration sequence autocorrelation peak; performing integrity verification on the multidimensional feature tensor through a convolutional neural network and outputting the integrity verification result; and performing smooth reconstruction on the local trajectory segments based on the integrity verification result to determine the final industrial robot motion trajectory representation.

[0082] In this embodiment, the input is a sequence of trajectory points corresponding to the obstacle avoidance optimal path. The trajectory points are arranged chronologically, and the sampling interval is 0.004 seconds (0.004 seconds is determined by the controller interpolation cycle to ensure that curvature calculation, velocity change detection, correlation calculation, and smooth resampling use the same time base). Each trajectory point includes three-dimensional position, three-dimensional velocity, and three-dimensional acceleration. The three-dimensional velocity is obtained by dividing the position difference of two adjacent trajectory points by 0.004 seconds, and the three-dimensional acceleration is obtained by dividing the velocity difference of two adjacent trajectory points by 0.004 seconds. The velocity and acceleration of the first trajectory point are a weighted combination of the first two differences, with weights of 0.6 and 0.4 (0.6 and 0.4 are jointly determined by boundary difference noise suppression and continuity requirements to ensure that the first point does not have numerical spikes and does not introduce future point information).

[0083] For obstacle avoidance path optimization, first calculate the set of curvature extrema points and the set of velocity change points, and then calculate the spatial overlap density based on these. The calculation of curvature extrema points is based on the bending strength of the position sequence, generated point by point: for the second trajectory point to the second-to-last trajectory point, take the previous point, the current point, and the next point to form a three-point neighborhood, and calculate the displacement vector of the previous segment and the displacement vector of the next segment respectively; each component of the displacement vector is obtained by subtracting the corresponding coordinates, and the displacement length is obtained by taking the square root of the sum of the squares of each component; divide each component of the displacement vector by the displacement length to obtain the unit direction vector; then calculate the dot product of the unit direction vector of the previous segment and the unit direction vector of the next segment, the dot product is obtained by multiplying the corresponding components and summing them one by one, and convert the dot product into the turning amplitude; normalize the turning amplitude according to the neighborhood average step size to obtain the bending strength, the neighborhood average step size is taken as the arithmetic mean of the displacement length of the previous segment and the displacement length of the next segment (the arithmetic mean is determined by the step size jitter suppression requirement to avoid abnormal amplification of bending strength by unilateral displacement). After the bending intensity sequence is generated, extreme value detection is performed. The detection window length is 5 trajectory points (5 is determined by the statistics of the shortest continuous width of the trajectory turning event sampled at 0.004 seconds. The statistical basis is the upper bound of the number of sampling points converted from the half-width of the turning peak in the historical trajectory, which is 2 and then covered by the center point to form 5). When the bending intensity of the center point of the window is strictly greater than that of the other 4 points in the window, it is marked as a curvature maximum point. When the bending intensity of the center point of the window is strictly less than that of the other 4 points in the window, it is marked as a curvature minimum point. The two types of points are merged to form a curvature extreme point set. To eliminate noise triggering in straight sections, an effective lower limit for bending strength is set. The maximum bending strength of a straight section is multiplied by 3 (the multiplier of 3 is determined by the interval between the upper limit of noise in straight sections and the lower limit of bending strength in turning sections, so that the false triggering rate of straight sections is reduced to below 0.10 and the detection rate of turning sections is maintained above 0.90; the straight section is determined by the interval where the acceleration amplitude is continuously lower than the steady-state threshold and the duration is not less than 32 trajectory points. The steady-state threshold is the upper limit of servo steady-state acceleration jitter multiplied by 1.2, and 1.2 is determined by the upper limit of jitter increment caused by load change). Extreme points where the bending strength is lower than this lower limit are not included in the set.

[0084] The calculation of velocity abrupt change points is generated point-by-point based on the velocity amplitude change. For each trajectory point, the velocity amplitude is calculated, obtained by summing the squares of the three-dimensional velocity components and then taking the square root. For the second to the last trajectory point, the velocity change amplitude is calculated, taken as the absolute value of the current velocity amplitude minus the velocity amplitude of the previous point. The velocity abrupt change threshold is the maximum velocity change amplitude of the uniform velocity segment multiplied by 1.1 (1.1 is determined by superimposing the upper bound of velocity measurement noise and the upper bound of interpolation error, used to cover batch fluctuations and suppress false triggering by noise; the uniform velocity segment is determined by an interval where the acceleration amplitude is continuously lower than the steady-state threshold for a duration of no less than 32 trajectory points, and the steady-state threshold is determined in the same way as before). When the velocity change amplitude of a trajectory point is strictly greater than this threshold, the point is marked as a velocity abrupt change point, forming a set of velocity abrupt change points.

[0085] Spatial overlap density is used to measure the spatial strength of the correspondence between curvature extrema and velocity mutation points. For each curvature extrema, velocity mutation points within two trajectory points before and after the time index are selected as candidate sets (the two points are determined by the upper bound of the maximum offset of event alignment, which is obtained by statistically analyzing historical data to ensure that the time misalignment between the curvature extrema and velocity mutation does not exceed two sampling points). For each velocity mutation point in the candidate set, the spatial distance between it and the curvature extrema is calculated. The spatial distance is obtained by taking the square root of the sum of the squared differences of the three-dimensional coordinates. The spatial matching tolerance is 20 mm (20 mm is determined by adding the upper bound of the repeatability error, the upper bound of the perception error, and the margin of the end-shape boundary to ensure that the matching criterion covers the trajectory drift). When there are candidate points with a distance of no more than 20 mm, the one with the smallest distance is selected to establish a pair of matches and counted as 1 overlap; when there are no candidate points that meet the tolerance, it is counted as 0 overlaps. To avoid the same velocity mutation point being occupied repeatedly, a one-to-one matching constraint is adopted, and the matched velocity mutation point is removed from the subsequent candidate set. The spatial overlap density is obtained by summing up the number of times all curvature extrema coincide and dividing by the number of curvature extrema. When the number of curvature extrema is 0, the spatial overlap density is 0 to prevent false high density from appearing in segments without events.

[0086] Local trajectory segments are extracted based on spatial overlap density, and the range of segments requiring integrity verification is determined. First, the temporal index of each pair of matching points is taken as the key index. Local trajectory segments are then extracted from the obstacle avoidance optimization path based on the key index. For each local trajectory segment, 64 trajectory points are taken forward and 64 trajectory points are taken backward from the key index, forming a local segment with a length of 129 trajectory points (64 points are determined by the local structure coverage and the receptive field requirements of the subsequent convolutional network, based on the principle that the longest continuous width of a turning event from entry to exit does not exceed 128 sampling points; 129 points are obtained by adding the center point to the previous and next 64 points, ensuring that the key point is located at the center of the local segment). When the key index is close to the start or end point of the trajectory, resulting in fewer than 64 points, boundary mirroring is used for padding. The length of the mirroring padding is the number of missing points. Mirroring padding is completed by copying existing points within the boundary in reverse order to the missing side, ensuring that the local segment length is fixed at 129 points. Spatial overlap density is used as a triggering condition, and the triggering threshold is set to 0.30 (0.30 is determined by historical trajectory statistics: the spatial overlap density is compared with the subsequent geometric structure destruction labels, the scanning threshold is increased by 0.05 from 0.10 to 0.60, and the threshold with a false trigger rate of no more than 0.10 and a missed trigger rate of no more than 0.10 is selected to obtain 0.30). When the spatial overlap density is not less than 0.30, the local segment of the corresponding key index is included in the integrity verification set.

[0087] For each local trajectory segment included in the set, geometric self-similarity and the peak value of the acceleration sequence autocorrelation are calculated, and a multidimensional feature tensor is constructed accordingly. Geometric self-similarity is used to measure the morphological consistency of local segments at different scales, and the calculation is completed in the order of "arc length normalization, resampling, alignment, and error mapping". First, the arc length sequence is calculated for the local segment position sequence: the spatial distance between two adjacent position points within the local segment is calculated, and the spatial distance is obtained by summing the squares of the differences in the three-dimensional coordinates and then taking the square root; the adjacent distances are accumulated segment by segment to form the cumulative arc length, and the cumulative arc length is divided by the total arc length to obtain the normalized arc length. Then, two scale sequences are constructed: scale 1 takes the original 129-point sequence, and scale 2 takes 1 point out of every 2 points to form a sparse sequence (the choice of 2 is determined by the constraint that the local morphology still maintains the turning contour at this scale and the amount of computation is controlled, and the constraint is based on the fact that the key turning segment still covers no less than 20 points after sparsening). Resampling was performed on both scale 1 and scale 2 under the same arc length coordinates. The number of resampling points was 64 (64 was determined by the shape comparison resolution and the convolution input length; after comparing candidates 48, 64, and 80, 64 was the optimal score in both error variance and computation delay). The resampling method involved interpolating points at 64 fixed positions with normalized arc lengths from 0 to 1, using linear interpolation within adjacent arc length intervals. Rigid alignment was performed on the two 64-point sequences. The alignment process first calculated the centroid of each sequence, which was obtained by summing the coordinates of the 64 points and dividing by 64. The centroid of each sequence was then subtracted from the coordinates of each point in the sequence to complete the translation alignment. Next, the principal direction was calculated and rotation alignment was performed. The principal direction was determined by the direction of maximum dispersion of the point coordinate covariance. The covariance was calculated by multiplying and summing the points along the coordinate dimension after removing the mean. The direction of maximum dispersion was obtained through symmetric matrix decomposition and the direction of maximum variance contribution was selected. The two sequences were rotated until their principal directions were consistent to complete the alignment. After alignment, the distance error sequence of corresponding points is calculated. The error is obtained by summing the squares of the differences between the coordinates of corresponding points in the two sequences and then taking the square root. The 64 errors are squared point by point and summed, then divided by 64 to obtain the mean square error. The mean square error is mapped to a geometric self-similarity score. The mapping adopts an inverse linear mapping: the score is 1 when the mean square error is 0, the score is 0 when the mean square error reaches the upper limit of the error, and the score is 0 when it exceeds the upper limit of the error. The upper limit of the error is the maximum mean square error of the alignment of the local segment of the strongly coupled sample multiplied by 1.10 (1.10 is determined by the upper limit of the operating condition fluctuation. The strongly coupled samples are derived from the previous strongly coupled screening and have been manually checked for consistency. The multiplier of 1.10 ensures that the strongly coupled local segment still falls within the positive score interval under the field noise).

[0088] The autocorrelation peak value of the acceleration sequence is used to measure the repeatability of the dynamic rhythm in a local segment. The calculation is completed in the following order: amplitude sequence generation, mean removal, autocorrelation accumulation, and peak extraction. First, the acceleration amplitude is calculated for each trajectory point in the local segment. The acceleration amplitude is obtained by summing the squares of the three-dimensional acceleration components and then taking the square root, forming an amplitude sequence of length 129. The amplitude sequence is then mean-removed. The mean is obtained by summing the 129 amplitudes point by point and then dividing by 129. Mean removal is done by subtracting the mean from the amplitude at each point. The lag range is set to 1 to 32 trajectory points (32 is converted from the maximum period of a typical acceleration / deceleration rhythm at a sampling time of 0.004 seconds, and the statistical basis is that the distance between the main acceleration peaks in the historical trajectory does not exceed 32 sampling points). Autocorrelation accumulation is performed on each lag value: the mean-removed sequence is multiplied point by point by its lag version and accumulated. The number of accumulated terms is 129 minus the lag value. The accumulated result is then divided by the number of accumulated terms to obtain the normalized autocorrelation value. After traversing all lag values, the maximum positive autocorrelation value is taken as the autocorrelation peak value, and the corresponding lag position is recorded as the peak position feature; when the maximum positive autocorrelation value is less than 0, the autocorrelation peak value is set to 0 and the peak position feature is set to 0.

[0089] After obtaining the geometric self-similarity score and autocorrelation peak, a multidimensional feature tensor is constructed and fed into a convolutional neural network to perform integrity-preserving verification. The time dimension of the multidimensional feature tensor is fixed at 64 (64 is consistent with the aforementioned number of resampling points, ensuring that geometric and dynamic features are aligned by the same index), and the number of channels is fixed at 4 and arranged in a fixed order. The first channel is the bending intensity sequence, which is obtained by normalizing the bending intensity of local segments by arc length and then resampling to 64 points; the bending intensity calculation method is consistent with that used for curvature extremum detection. The second channel is the velocity change amplitude sequence, which is obtained by first calculating the velocity amplitude of local segments point by point, then obtaining the velocity change amplitude sequence by the absolute value of the difference between the velocity amplitudes of two adjacent points, and then resampling to 64 points after normalizing by arc length. The third channel is the acceleration amplitude sequence, which is obtained by obtaining 129 acceleration amplitudes using the aforementioned method and then resampling to 64 points after normalizing by arc length. The fourth channel is a global scalar channel. Geometric self-similarity scores and autocorrelation peaks are first combined into a single scalar at a predetermined ratio of 0.6 and 0.4 (these ratios are determined by the statistical weights of the contributions of geometric and dynamic damage to the integrity label; the statistical method involves calculating the sensitivity of each index to the misclassification rate on the labeled data and normalizing the weights). This scalar is then copied 64 times to form a sequence of length 64, which serves as the fourth channel. All four channels are standardized before entering the network. The standardized mean and standardized standard deviation are obtained from the modeling phase. The standardized mean is obtained by summing the training samples point by point and dividing by the number of samples. The standardized standard deviation is obtained by removing the mean, summing the squares point by point, dividing by the number of samples minus 1, and then taking the square root. The lower limit of the standardized standard deviation is 0.001 (0.001 is determined by numerical stability; values ​​below 0.001 will amplify small noise and create spikes in the convolution output). When the standard deviation of a channel is less than 0.001, 0.001 is used in the division.

[0090] The convolutional neural network performs integrity-preserving verification calculations on multidimensional feature tensors in a hierarchical order. The network employs three 1D convolutional layers, one global average pooling layer, and two fully connected classification heads. The first convolutional layer has a kernel length of 5 and 32 output channels (the kernel length of 5 is determined by the local event width coverage requirement, statistically based on the fact that the peak bending intensity does not exceed 5 sampling points in the full width at half maximum after resampling; the number of channels 32 is determined by the constraints of classification accuracy and inference latency, choosing 32 from candidates 16, 32, and 48 to ensure that the single-segment inference latency does not exceed 1 millisecond and achieves the highest accuracy). The second convolutional layer has a kernel length of 5 and 64 output channels (64 is determined by the complexity of the mid-layer feature combination, chosen to reduce the false positive rate to below 0.10 under the same latency budget). The third convolutional layer has a kernel length of 3 and 64 output channels (the kernel length of 3 is used to enhance the discrimination of edge abrupt changes, chosen to improve the detection rate of edge differences before and after local reconstruction). Each convolutional operation is performed within a sliding time window. Within the window, the input and convolution weights are multiplied item by item, summed, and then the bias is added to obtain the convolutional output. The convolution stride is 1 (determined by maintaining temporal resolution to avoid event position drift caused by downsampling). Symmetrical padding is used, with the padding width matching the kernel length to maintain an output length of 64 (the padding strategy is determined by the consistency requirement of input and output lengths to ensure that the size of subsequent pooling and the classification head remains fixed). Each convolutional output undergoes linear rectification, which sets negative values ​​to 0 to suppress invalid negative responses. Global average pooling accumulates each output channel along the 64-point time dimension and divides by 64 to obtain the channel summary vector. The first fully connected layer of the classification head maps the channel summary vector to a hidden vector of length 128, the length of which is determined by both the expressive power of the classification boundary and the computational cost. The second fully connected layer outputs two confidence scores, corresponding to "intact" and "destructive," respectively. The confidence scores are obtained through a normalization function transformation, ensuring that the sum of the two is 1. The integrity judgment threshold is set to 0.70 (0.70 is determined by threshold scanning: the threshold is traversed from 0.50 to 0.90 in increments of 0.05 on the labeled data, and the threshold with a false positive rate of no more than 0.10 and a false negative rate of no more than 0.10 is selected, and finally 0.70 is obtained). When the "integrity" confidence is not less than 0.70, the integrity verification is considered to have passed; otherwise, it is considered to have failed.

[0091] Based on the integrity verification results, smoothing reconstruction is performed on local trajectory segments to generate the final industrial robot motion trajectory representation. For verified local segments, a light smoothing is applied, which only affects the position sequence and keeps the endpoint positions unchanged. The smoothing window length is set to 5 trajectory points (5 is determined by the combined constraints of noise suppression and peak position preservation; when set to 5, the bending extreme value time index offset does not exceed 1 sampling point). For the midpoint of the local segment, the average value of each point within the window is used to replace the original position point. The average value is calculated by summing the three-dimensional coordinates and then dividing by 5. For the two points at each end of the local segment, mirror completion is used to participate in the average value calculation, with the mirror completion length set to 2 to match the window radius. After smoothing, the velocity and acceleration are recalculated with a difference of 0.004 seconds. The jerk is then obtained from the acceleration difference, and the jerk amplitude is calculated. The upper limit of the jerk amplitude is the maximum allowable jerk amplitude of the controller multiplied by 0.90 (0.90 is determined by the safety margin to cover the instantaneous overshoot caused by load changes and friction changes). When the jerk amplitude at any point in a local segment exceeds the upper limit, the refactoring process is used instead of the light smoothing to avoid directly outputting the non-executable segment.

[0092] For local segments that fail verification, constraint spline reconstruction is performed. The reconstruction uses cubic splines and simultaneously constrains endpoint positions, endpoint tangentials, accelerometer amplitude, and obstacle safety distance. Endpoint position constraints ensure the first and last points of the local segment remain unchanged. Endpoint tangential constraints ensure the first and last tangentials are consistent with the original local segment; the tangential is obtained by normalizing the position differences between the endpoint and its adjacent points, with the normalization method being dividing each component of the difference vector by the difference length. Smoothing constraints use the upper limit of accelerometer amplitude, with the value determined in the same way as for light smoothing. The safety distance constraint is 20 mm (20 mm is determined by the obstacle safety envelope expansion distance to ensure local reconstruction does not disrupt obstacle avoidance gaps). Distance calculation is performed based on the closest distance from the trajectory point to the safety envelope, and the closest distance calculation method is consistent with near-obstacle recognition. The spline control points are solved using iterative updates. The update objective is to simultaneously satisfy the above constraints and minimize the sum of squares of position changes. The sum of squares of position changes is obtained by calculating the difference between the reconstructed position and the original position point by point and then summing the squares. The iteration stopping condition is that the decrease in constraint violation between two adjacent iterations does not exceed 0.00001 (0.00001 is determined by both numerical accuracy and execution error sensitivity; improvements below this value cannot be distinguished from sampling noise). The maximum number of iterations is 50 (50 is determined by the planning cycle delay budget and the statistics of the slowest convergence sample to ensure that real-time performance is still met in the worst case). After reconstruction, 129 local segments are generated by resampling at 0.004 seconds. The safety distance of each resampled point is checked to be no less than 20 mm and the jerk amplitude does not exceed the upper limit. After all points meet the requirements, the reconstructed local segment is used as the replacement segment.

[0093] When multiple local segments overlap at the time index, weighted fusion is used to eliminate splicing abrupt changes. For each time index position within the overlapping area, the distance to the center index of the local segment is calculated. When the distance is 0, the center weight is 1; when the distance reaches 64, the boundary weight is 0; and the weight of the intermediate positions is calculated linearly decreasing. This weight is used to perform a weighted summation of the reconstructed segment position and the original obstacle avoidance preferred path position to obtain the fused position. Based on the fused position, the velocity and acceleration are re-differentiated to ensure that the fused sequence still maintains time alignment. After all local segments are replaced and fused, the output full trajectory sequence serves as the final representation of the industrial robot's motion trajectory.

[0094] S8 includes: extracting the geometric symmetry axis of the final industrial robot motion trajectory representation and calculating the axis alignment at the moment of velocity vector flipping of the original trajectory data; extracting the curvature zero-crossing frequency of the trajectory representation and calculating the frequency matching degree with the velocity zero-crossing frequency of the original trajectory data; calculating the consistency score based on the axis alignment and frequency matching degree; if the consistency score is lower than a preset threshold, generating an error feedback signal to update the joint principal component axis weights in the initial representation vector, and re-entering the dimensionality reduction process based on the updated joint principal component axis weights.

[0095] In this embodiment, the final industrial robot motion trajectory representation is first aligned with the original trajectory data by time index before entering the consistency evaluation. The sampling interval is 0.004 seconds (0.004 seconds is determined by reading the controller interpolation cycle, which is a controller operating parameter; this value ensures that the differential calculation of trajectory points, velocity, and acceleration, as well as event detection, are all on the same time base). The alignment verification includes two items: the difference in the starting timestamp does not exceed one sampling interval, and the difference in the length of the two sequences does not exceed one sampling point (the tolerance of one sampling point is determined by the statistics of communication jitter and log landing point error; a tolerance exceeding one sampling point will introduce a mismatch in the speed flip time). If the verification is not met, resampling is performed. Resampling uses the timestamp as an index and linear interpolation is used to fill in the missing positions. The interpolation span does not exceed two sampling intervals (two sampling intervals are determined by the statistical upper limit of the maximum power transmission jitter on site; a span exceeding two will significantly amplify the error in the speed direction change and affect the flip judgment).

[0096] The geometric symmetry axis is extracted from the 3D position sequence represented by the final motion trajectory. The extraction process is executed according to the following steps: centering, candidate axis generation, mirror coincidence evaluation, and optimal axis determination. In the centering stage, the centroid is first calculated. The three coordinate components of the centroid are obtained by summing the corresponding components of all position points and dividing by the number of points. Then, the centroid is subtracted from each position point to obtain the centered point set. Centering is used to eliminate the interference of translation on the evaluation of the symmetry axis direction. In the candidate axis generation stage, the covariance matrix of the centered point set is calculated. The covariance matrix is ​​obtained by multiplying each coordinate component pairwise after removing the mean and summing the results over all points. Then, it is divided by the number of points and subtracted by 1 to complete the normalization (subtracting 1 from the number of points in the denominator is used to reduce the systematic bias of underestimating the variance of finite samples). The covariance matrix is ​​symmetrically decomposed to obtain several principal directions, which are then sorted in descending order of dispersion. The top 8 principal directions are selected as candidate axes (these 8 directions are jointly determined by the search resolution and the evaluation delay budget: fewer than 8 candidate axes will result in a significant increase in the deviation of the symmetry axis direction; more than 8 candidate axes will significantly increase the number of nearest point searches in the mirror coincidence evaluation and exceed the planning cycle budget). In the mirror coincidence evaluation stage, point-by-point mirroring is performed for each candidate axis, and the average mirror distance is calculated: for any centralized position point, the perpendicular foot point to the candidate axis is first calculated. The perpendicular foot point is obtained by multiplying the projection length of the position point in the axis direction by the axis direction component and adding it to the reference point on the axis; then, the vertical component of the point relative to the perpendicular foot point is inverted to form the mirror point; the nearest point to the mirror point is searched in the centralized point set, and the distance between the two is calculated. The distance is obtained by taking the square root of the sum of the squared differences of the three-dimensional coordinates. Closest point search employs a spatial partition index to control computational load. The partition side length is set to 20 mm (this 20 mm is determined by the average spacing between trajectory points and the real-time search overhead: too small a partition side length leads to an excessive number of partitions and increased index maintenance costs; too large a partition side length leads to an excessive number of candidate points and increased closest point traversal times; 20 mm ensures a stable number of candidate points and meets the search latency budget under the current trajectory point density). The average mirror distance is obtained by summing the closest distances corresponding to all mirror points and dividing by the number of points. A smaller average mirror distance indicates a higher degree of mirror overlap. A reverse mapping is performed on the average mirror distance to obtain a symmetry score, where the score is 1 when the average mirror distance is 0, and 0 when the average mirror distance reaches the upper limit. The upper limit is taken as the average point spacing of the trajectory segment multiplied by 5 (the average point spacing is obtained by summing the distances of adjacent points and dividing by the number of adjacent segments; the multiplier of 5 is determined by the mirror matching tolerance coverage requirements, and with a multiplier of 5, slight asymmetry still shows score differences and is not prematurely suppressed to 0). After calculating the symmetry score for all candidate axes, the candidate axis with the highest score is selected as the geometric symmetry axis, and its direction vector is normalized. Normalization is accomplished by dividing each component of the direction vector by its length. The length is obtained by summing the squares of the components and taking the square root. Normalization is used to ensure that the scale of subsequent angle calculations is consistent.

[0097] The velocity vector flipping moment is detected from the original trajectory data. The detection process is executed according to the following steps: velocity vector generation, noise suppression, flipping detection, and persistence confirmation. In the velocity vector generation stage, if the original data contains a velocity vector, it is directly used; if the original data only contains position, the velocity vector is obtained by difference: the current position minus the previous position yields the displacement vector, and each component of the displacement vector is divided by 0.004 seconds to obtain the velocity vector. The boundary point velocity is obtained by a weighted combination of the differences of two adjacent segments, with weights of 0.6 and 0.4 (0.6 and 0.4 are determined by boundary difference noise suppression: 0.6 increases the contribution of the nearest segment to maintain true local changes, and 0.4 suppresses the initial peak caused by single-segment anomalies). In the noise suppression stage, a moving average filter is applied to each component of the velocity vector, with a window length of 5 sampling points (5 is determined by statistics of the velocity jitter frequency and the duration width of the flipping event: a window length less than 5 results in jitter residue triggering false judgments; a window length greater than 5 increases the offset of the flipping moment and affects axis alignment calculations). In the flipping discrimination stage, the angle between the current velocity direction and the previous velocity direction is calculated for each sampling point. The angle is obtained by converting the dot product and length relationship: the dot product is obtained by multiplying the corresponding components and then summing them one by one, and the length is obtained by summing the squares of the components and then taking the square root. When the angle is greater than 150 degrees and the velocity amplitude at both moments is not lower than the lower velocity limit, it is marked as a flipping candidate (150 degrees is determined by the lower bound of the angle distribution of real reverse events in historical data, and taking 150 degrees can exclude the direction change caused by ordinary turning; the lower velocity limit is the average velocity amplitude of the uniform speed segment multiplied by 0.2, and 0.2 is determined by the upper bound of the proportion of low speed noise. When it is below this threshold, the velocity direction is easily dominated by noise and has no discrimination significance; the uniform speed segment is determined by the interval where the acceleration amplitude is continuously lower than the steady-state threshold and the duration is not less than 32 sampling points, and 32 is determined by the statistics of the shortest duration of the stable segment; the steady-state threshold is the upper bound of the servo steady-state acceleration jitter multiplied by 1.2, and 1.2 is determined by the upper bound of the jitter increment caused by load change). The continuous confirmation phase requires that a flip candidate continuously satisfy at least two sampling points before it is confirmed as a flip moment (the two sampling points are statistically determined by the possibility of single-point false flips caused by interpolation jitter; the continuity condition can significantly reduce false detections). The set of velocity vector flip moments is obtained after confirmation.

[0098] Axis alignment is used to measure the consistency of the direction between the geometric symmetry axis and the velocity reversal event. The calculation is performed by accumulating each reversal event and averaging the results. For each reversal moment, the velocity direction vectors of the sampling point before and after the reversal are taken and normalized. The normalization method is to divide each component by the vector length. The principal reversal direction is obtained by adding the direction vectors before and after the reversal, and the sum is normalized. When the length of the sum is less than 0.1, the difference vector between the two direction vectors is used as the principal reversal direction and normalized (0.1 is determined by the lower limit of numerical stability of normalization; a sum length less than 0.1 will amplify numerical errors and lead to instability in the principal direction). The angle between the main flip direction and the geometric symmetry axis is then calculated and mapped to an alignment score: 1 for an angle no greater than 15 degrees, 0.7 for an angle greater than 15 degrees but no greater than 45 degrees, 0.4 for an angle greater than 45 degrees but no greater than 90 degrees, and 0 for an angle greater than 90 degrees (15 degrees, 45 degrees, and 90 degrees are determined by statistically analyzing the angle distribution boundaries between consistent and inconsistent samples in historical trajectories, selecting the latter two types of samples with the least overlap). The alignment scores of all flip events are summed and divided by the number of flip events to obtain the axis alignment degree; when the number of flip events is 0, the axis alignment degree is 0.

[0099] Frequency matching degree is used to measure the consistency between the curvature zero-crossing frequency of the final trajectory and the velocity zero-crossing frequency of the original trajectory. The calculation employs window-based statistical analysis and alignment matching. The curvature zero-crossing frequency is extracted from the final trajectory: first, the curvature sign and curvature intensity are calculated point-by-point. For the second to the second-to-last point, a three-point neighborhood is taken, and the unit direction of the preceding and following segments is calculated, with the curvature intensity obtained from the turning amplitude. The curvature intensity is normalized according to the average step size of the neighborhood. The curvature sign is determined according to a fixed reference direction to ensure consistent signs across the entire trajectory. The curvature signs are traversed along time; a zero-crossing is counted when two adjacent points have opposite signs and both points have curvature intensities higher than the lower limit of curvature noise. The lower limit of curvature noise is taken as the maximum bending intensity of the straight segment multiplied by 3 (the multiplier of 3 is statistically determined by the interval between the upper limit of straight noise and the lower limit of the actual steering bending intensity; after taking 3, the false trigger rate of the straight segment is not higher than 0.1 and the detection rate of the steering segment is not lower than 0.9; the straight segment is determined by the interval where the acceleration amplitude is continuously lower than the steady-state threshold and the duration is not less than 32 sampling points; the steady-state threshold is taken as the upper limit of servo steady-state acceleration jitter multiplied by 1.2, and the determination of 1.2 is based on the same basis as before). The zero-crossing frequency of velocity is extracted from the original trajectory: first, calculate the velocity amplitude sequence; calculate the mean velocity amplitude in each statistical time window and subtract the mean from the velocity amplitude at each point to obtain the velocity deviation sequence; traverse the signs of the velocity deviations, and when two adjacent points have opposite signs and the deviation amplitudes of both points are higher than the lower limit of velocity noise, it is counted as one zero-crossing. The lower limit for velocity noise is taken as the maximum velocity deviation of the uniform velocity segment multiplied by 1.1 (1.1 is determined by superimposing the upper limit of velocity measurement noise and the upper limit of interpolation error, used to cover batch fluctuations and suppress noise zero-crossing; the uniform velocity segment is determined by the interval where the acceleration amplitude is continuously lower than the steady-state threshold and the duration is not less than 32 sampling points). The same time window is used for both types of frequency statistics, with a time window length of 64 sampling points and a step size of 16 sampling points (64 and 16 are jointly determined by the frequency count variance and response delay: 64 significantly reduces the zero-crossing count variance, and 16 ensures sufficient overlap between adjacent windows, forming a smooth sequence while keeping the computational load under control). The number of curvature zero-crossings and the number of velocity zero-crossings are counted for each time window, and each is divided by the coverage duration of that window to obtain the frequency. The coverage duration is obtained by multiplying the number of adjacent sampling segments, 63, by 0.004 seconds. For each time window, the absolute value of the difference between two frequencies is calculated and then normalized to the upper limit of acceptable difference, and inversely mapped to the matching score for that window. The score is set so that a difference of 0 results in a score of 1, a difference reaching the upper limit results in a score of 0, and a difference exceeding the upper limit results in a score of 0. The upper limit of acceptable difference is the maximum value of the two frequency differences in the consistent samples multiplied by 1.1 (1.1 is determined by the upper bound of the operating condition fluctuation; consistent samples are derived from historical trajectory records that have passed consistency judgment and have been manually sampled and confirmed). The frequency matching degree is obtained by summing the matching scores of all time windows and dividing by the number of time windows.

[0100] The consistency score is obtained by weighted summation of axis alignment and frequency matching, with weights of 0.5 and 0.5 (0.5 and 0.5 are determined by balancing the contribution of validation data: the weight combinations are iterated over on the labeled data, constraining the false positive rate to be no higher than 0.1 and the false negative rate to be no higher than 0.1, and the combination with the highest accuracy is selected to obtain equal weights). The consistency threshold is 0.7 (0.7 is determined by threshold scanning: the threshold is iterated over on the labeled data from 0.5 to 0.9 in increments of 0.05, filtering for thresholds with a false positive rate no higher than 0.1 and a false negative rate no higher than 0.1, and the one with the highest accuracy is selected to obtain 0.7). When the consistency score is lower than 0.7, an error feedback signal is generated and the joint principal component axis weights in the initial representation vector are updated.

[0101] The error feedback signal is the difference between the consistency threshold of 0.7 and the current consistency score; a positive difference indicates insufficient consistency. The joint principal component axis weights correspond to the contribution coefficients of each joint principal axis within the reduced-dimensional subspace. The weights are non-negative and sum to 1. Weight updates are performed using sensitivity estimation, weight correction, normalization, and amplitude limiting. Sensitivity estimation employs perturbation evaluation: for each joint principal axis, its weight is increased by 0.01, while the remaining weights are proportionally reduced to maintain a sum of 1 (0.01 is determined by numerical stability and evaluation overhead: increments less than 0.01 are easily drowned out by computational noise, while increments greater than 0.01 cause weight structure jumps and amplify nonlinear effects). Subsequently, a dimensionality reduction mapping, trajectory reconstruction, and consistency score recalculation are performed with the updated weights. The sensitivity is the difference between the old and new consistency scores divided by 0.01. The weight correction is calculated by multiplying the error feedback signal by the sensitivity, and then multiplying by the update step size of 0.1 (0.1 is determined by both the convergence speed and weight oscillation suppression: a step size greater than 0.1 will cause the weight to oscillate more between adjacent updates, while a step size less than 0.1 will cause the update cycle to increase and exceed the planned cycle budget). After updating the weights for each joint principal axis according to its correction amount, non-negativity constraints are applied, and weights with negative values ​​are directly set to 0. Then, all weights are summed to obtain the total, and each weight is divided by this total to normalize it, restoring the total weight sum to 1. To avoid numerical amplification caused by an excessively small total sum, if the total sum is less than 0.001, the weights are reset to a uniform distribution (0.001 is determined by the lower limit of numerical stability; values ​​below this value will cause normalization amplification errors). After the weights are updated, a dimensionality reduction process is performed according to the updated joint principal component axis weights to generate a new dimensionality-reduced joint feature representation, trajectory reconstruction, and consistency evaluation results. The maximum number of updates is 5 (5 is determined by the planning cycle delay budget and the cost of recalculating a single update; even in the worst case, the real-time constraint is still satisfied when the upper limit is 5). When the consistency score reaches or exceeds 0.7, the update stops and the current final industrial robot motion trajectory representation is output.

[0102] Although embodiments of the invention have been shown and described, it will be understood by those skilled in the art that various changes, modifications, substitutions and alterations can be made to these embodiments without departing from the principles and spirit of the invention, the scope of which is defined by the appended claims and their equivalents.

Claims

1. A method for obstacle avoidance in the motion trajectory of an industrial robot based on deep learning, characterized in that, include: S1. Extract position sequence, velocity sequence, and acceleration sequence from the trajectory data of the industrial robot's motion trajectory, generate geometric feature vector and dynamic feature vector respectively through a convolutional neural network, and fuse the two to obtain an initial representation vector; S2. Calculate the synchronization between local curvature change and acceleration peak, spatial overlap density between curvature extreme points and velocity change points, and directional consistency between trajectory torsion gradient and acceleration rate of change vector for the initial representation vector. If the weighted average of the three indicators exceeds the preset threshold, it is determined that there is a strong coupling relationship between geometry and dynamic characteristics. S3. Based on the initial representation vector determined to be strongly coupled, principal component analysis is used to extract the geometric-dynamic coupling principal direction that retains the maximum variance, the joint principal component axis in the cross-modal information set, and the subspace that minimizes the joint reconstruction error of geometric and dynamic distortion after projection, so as to obtain the dimension-reduced joint feature representation. S4. Reconstruct the sequence of industrial robot trajectory sample points from the dimensionality-reduced joint feature representation, input it into the long short-term memory network to capture the temporal pattern corresponding to the matching degree of curvature zero-crossing frequency and velocity zero-crossing frequency, as well as the matching trend of the distribution of bending segment length and the statistical distribution of acceleration amplitude, and obtain the trajectory dynamics pattern vector.

2. The obstacle avoidance method for industrial robot motion trajectory based on deep learning according to claim 1, characterized in that: S1 includes: Acquire trajectory data during the operation of an industrial robot, and parse the trajectory data to obtain position sequence, velocity sequence and acceleration sequence; The position sequence is input into a spatial feature extraction convolutional neural network to generate a geometric feature vector, and the velocity sequence and acceleration sequence are input into a dynamic feature extraction convolutional neural network to generate a dynamic feature vector. The geometric shape feature vector and the dynamic feature vector are dimensionally mapped and feature concatenated to obtain an initial representation vector containing complete spatiotemporal information.

3. The obstacle avoidance method for industrial robot motion trajectory based on deep learning according to claim 1, characterized in that: S2 includes: The geometric and dynamic feature components are obtained by analyzing the initial representation vector. The local curvature change sequence and acceleration peak sequence are then extracted from the geometric and dynamic feature components. Calculate the degree of synchronization between the local curvature change sequence and the acceleration peak sequence, and determine the spatial overlap density between curvature extreme points and velocity change points based on the degree of synchronization. Based on the spatial overlap density, the directional consistency index is calculated. If the comprehensive coupling coefficient generated by the synchronization degree value, spatial overlap density and directional consistency index meets the conditions, it is determined that there is a strong coupling relationship between the geometry and dynamic characteristics.

4. The obstacle avoidance method for industrial robot motion trajectory based on deep learning according to claim 1, characterized in that: S3 includes: Obtain the initial representation vectors selected by the strong coupling decision, and use the initial representation vectors to construct the joint covariance matrix; The joint covariance matrix is ​​decomposed to extract the geometric dynamic coupling principal directions that retain the maximum variance, and the joint principal component axes in the cross-modal information set are identified. An initial subspace is constructed based on the geometric dynamic coupling principal direction and the joint principal component axis, and the target subspace is obtained by minimizing the joint reconstruction error function. The initial representation vector is mapped to the target subspace to obtain the dimensionality-reduced joint feature representation.

5. The obstacle avoidance method for industrial robot motion trajectory based on deep learning according to claim 1, characterized in that: S4 includes: The joint feature representation after dimensionality reduction is obtained, and the sequence of industrial robot trajectory sample points is obtained by inverse projection reconstruction. The frequency matching degree between curvature and velocity zero-crossing frequency and the distribution trend of bending segment length and acceleration amplitude are calculated based on the sequence of industrial robot trajectory sample points. The frequency matching degree and the distribution matching trend are combined to generate a time-series feature sequence, and the time-series feature sequence is input into a long short-term memory network; Obtain the hidden layer state vector output by the Long Short-Term Memory network, and map the hidden layer state vector to obtain the trajectory dynamics pattern vector.

6. The obstacle avoidance method for industrial robot motion trajectory based on deep learning according to claim 1, characterized in that, It also includes S5, which calculates the correspondence score between the path fractal dimension and the dynamic energy dissipation rate for the trajectory dynamics mode vector and the geometric shape feature vector. If the correspondence score is lower than a preset threshold, the positions of key control points in the industrial robot trajectory sample point sequence are adjusted iteratively through gradient descent until the adjusted sequence meets the predefined kinematic continuity constraints, specifically including: Obtain the trajectory dynamics mode vector and geometric shape feature vector, and calculate the path fractal dimension and dynamic energy dissipation rate respectively; The mutual information value between the path fractal dimension and the dynamic energy dissipation rate is used as the corresponding score; If the correspondence score is lower than the preset threshold, the control point position is updated along the gradient vector that maximizes the correspondence score in order to reconstruct the corrected industrial robot trajectory sample point sequence. The third-order differential of the corrected industrial robot trajectory sample point sequence is used to obtain the jerk sequence. If the jerk sequence meets the smoothness limit value, it is determined that the adjusted sequence meets the predefined kinematic continuity constraint.

7. The obstacle avoidance method for industrial robot motion trajectory based on deep learning according to claim 6, characterized in that, This also includes S6, which calculates the interference probability caused by the synchronization between local curvature changes and peak acceleration in the potential obstacle area based on the adjusted industrial robot trajectory sequence that satisfies kinematic continuity constraints, combined with the obstacle model in the industrial robot's workspace. If the interference probability is higher than a preset threshold, a dynamic programming algorithm is used to search for candidate node sequences that avoid the interference probability area, thus obtaining the optimal obstacle avoidance path. Specifically, this includes: Obtain the trajectory sequence and obstacle model of the industrial robot that satisfy the kinematic continuity constraints, and identify the potential obstacle region located near the obstacle model in the trajectory sequence; The degree of synchronization is calculated based on the local curvature changes and peak acceleration within the potential obstacle area to determine the probability of interference. If the interference probability is higher than a preset threshold, a dynamic programming algorithm is used to search for the optimal state transition path passing through candidate nodes outside the potential obstacle area to obtain the preferred obstacle avoidance path.

8. The obstacle avoidance method for industrial robot motion trajectory based on deep learning according to claim 7, characterized in that, This also includes S7, extracting key node locations with high spatial overlap density between curvature extrema and velocity mutation points from the obstacle avoidance optimization path, verifying the integrity of geometric self-similarity and acceleration sequence autocorrelation peaks through convolutional neural networks, and determining the final industrial robot motion trajectory representation, specifically including: Calculate the spatial overlap density between curvature extrema and velocity mutation points in the obstacle avoidance optimal path, and extract local trajectory segments based on the spatial overlap density; Obtain the geometric self-similarity and acceleration sequence autocorrelation peak of local trajectory segments, and construct a multidimensional feature tensor based on the geometric self-similarity and acceleration sequence autocorrelation peak; The integrity verification of the multidimensional feature tensor is performed using a convolutional neural network, and the integrity verification result is output. Based on the integrity verification results, a smooth reconstruction is performed on the local trajectory segments to determine the final representation of the industrial robot's motion trajectory.

9. The obstacle avoidance method for industrial robot motion trajectory based on deep learning according to claim 8, characterized in that, It also includes S8, calculating the consistency score between the final industrial robot motion trajectory representation and the original high-dimensional trajectory data of the industrial robot, specifically the alignment degree of the geometric symmetry axis and the velocity vector at the moment of reversal, as well as the matching degree of the curvature zero-crossing frequency and the velocity zero-crossing frequency. If the consistency score is lower than a preset threshold, the joint principal component axis weights in the initial representation vector are updated through error feedback, and the process re-enters the dimensionality reduction stage, specifically including: The geometric symmetry axis of the final industrial robot motion trajectory representation is extracted, and the velocity vector flipping moment of the original trajectory data is calculated to determine the axis alignment. The frequency matching degree is calculated by extracting the curvature zero-crossing frequency of the trajectory representation and the velocity zero-crossing frequency of the original trajectory data.

10. The obstacle avoidance method for industrial robot motion trajectory based on deep learning according to claim 9, characterized in that: S8 further includes: A consistency score is calculated based on the degree of axis alignment and frequency matching. If the consistency score is lower than the preset threshold, an error feedback signal is generated to update the joint principal component axis weights in the initial representation vector, and the vector re-enters the dimensionality reduction process based on the updated joint principal component axis weights.

Images