A control method and system for a robot arm

Abstract

The present application relates to the technical field of mechanical arm control, and particularly relates to a mechanical arm control method and system. The method comprises the following steps: performing acoustic wave monitoring by deploying a micro acoustic wave sensor at a mechanical arm joint, obtaining a friction acoustic wave frequency domain structure, and then deducing a radial excursion index of the mechanical arm joint, and performing spatial interpolation processing on the stress drift, and outputting stress drift regression data; by analyzing the regression data, quantifying the accuracy deviation of the joint, and performing deviation adaptive training correction based on the accuracy deviation data; finally, designing a deviation correction architecture and sending it to the terminal to execute mechanical arm control, realizing accuracy optimization and automatic correction of the mechanical arm. The present application optimizes the control technology of the mechanical arm, making the control technology of the mechanical arm more accurate.

Description

Technical Field

[0001] This invention relates to the field of robotic arm control technology, and in particular to a robotic arm control method and system. Background Technology

[0002] Robotic arms are increasingly widely used in modern manufacturing, logistics, and medical fields. As a highly efficient automation tool, the advantages of robotic arms in high precision, high speed, and high load capacity have led to their widespread application in many complex tasks. However, during long-term operation, especially under high load and high speed conditions, robotic arms are susceptible to wear, aging, and environmental changes, leading to friction between joints, changes in force, and deviations in motion accuracy, thus affecting the overall performance and task execution effectiveness of the robotic arm. To ensure the stability and efficiency of robotic arms, real-time monitoring of their operating status and accurate assessment and correction of deviations have become urgent problems to be solved. However, traditional robotic arm control methods suffer from inaccurate analysis of joint motion accuracy loss, resulting in large control correction errors. Summary of the Invention

[0003] Therefore, it is necessary to provide a control method and system for a robotic arm to solve at least one of the above-mentioned technical problems.

[0004] To achieve the above objective, a control method for a robotic arm is provided, the method comprising the following steps:

[0005] Step S1: Deploy miniature acoustic sensors at the joints of the robotic arm and monitor the acoustic waves during the joint operation to obtain joint operation monitoring acoustic waves; identify and label the joint operation monitoring acoustic waves with friction acoustic waves to output the frequency domain structure of friction acoustic waves.

[0006] Step S2: Based on the frequency domain structure simulation of the friction sound wave, derive the radial displacement index of the robotic arm joint during operation, and then perform spatial interpolation processing of the force drift of the robotic arm joint to output force drift regression data; quantify the running accuracy deviation of the joint based on the force drift regression data to obtain accuracy deviation quantification data.

[0007] Step S3: Based on the precision deviation quantification data, perform adaptive training and correction of deviations during the operation of the robotic arm to obtain adaptive correction data. Then, design the deviation correction architecture to build the deviation correction architecture and send the deviation correction architecture to the terminal to execute the control of the robotic arm.

[0008] The present invention also provides a control system for a robotic arm, for executing the control method for the robotic arm described above, the control system comprising:

[0009] The friction acoustic wave recognition module is used to deploy miniature acoustic wave sensors at the joints of a robotic arm and monitor the acoustic waves during the operation of the robotic arm joints to obtain the joint operation monitoring acoustic waves; the joint operation monitoring acoustic waves are marked with friction acoustic wave recognition to output the frequency domain structure of friction acoustic waves;

[0010] The precision deviation quantization module is used to derive the radial displacement index of the robotic arm joint during operation based on the frequency domain structure simulation of the friction acoustic wave, and then perform spatial interpolation processing of the force drift of the robotic arm joint to output force drift regression data; based on the force drift regression data, the running precision deviation of the joint is quantified to obtain precision deviation quantization data.

[0011] The correction architecture design module is used to perform adaptive training and correction of deviations during the operation of the robotic arm based on the precision deviation quantification data, obtain deviation adaptive correction data, and then design the deviation correction architecture to build the deviation correction architecture. The deviation correction architecture is then sent to the terminal to execute the control of the robotic arm.

[0012] The beneficial effects of this invention lie in the ability to monitor the working status of the robotic arm joints in real time, particularly friction and load conditions, by deploying miniature acoustic sensors at the joints. The application of acoustic sensors can capture minute frictional sound waves, and by identifying and labeling these waves, the frictional characteristics of the robotic arm joints can be accurately analyzed. This not only helps to provide early warning of potential problems but also provides a scientific basis for subsequent accuracy correction, thereby ensuring the efficient operation and long-term stability of the robotic arm. By converting the frequency domain structure of the frictional sound waves into a radial runout index, the dynamic changes of the robotic arm joints can be simulated and derived more accurately. This process not only helps to analyze the force on the joints under different loads but also uses spatial interpolation technology for force drift to output regression data, quantifying the operational accuracy deviation of the robotic arm joints. This allows for a clear understanding of the deviations during operation and provides quantitative data support for subsequent accuracy correction. Adaptive correction based on the quantified accuracy deviation data allows for real-time adjustment of the robotic arm's motion control, enabling the robotic arm to self-correct deviations during operation. This adaptive training and correction process effectively reduces human intervention and improves the automation and intelligence level of the system. The design of the deviation correction architecture makes the correction process more systematic, allowing for real-time execution and adjustment at the terminal. This ensures the robotic arm maintains optimal working accuracy and stability in various complex environments, significantly improving its efficiency and precision. Therefore, this invention optimizes a traditional robotic arm control method, addressing the problem of inaccurate joint movement accuracy loss analysis leading to large control correction errors. It improves the accuracy of joint movement accuracy loss analysis and reduces robotic arm control correction errors. Attached Figure Description

[0013] Figure 1 A flowchart illustrating the steps of a control method for a robotic arm;

[0014] Figure 2 for Figure 1 A detailed flowchart illustrating the implementation steps of step S2.

[0015] Figure 3 for Figure 1 A detailed flowchart illustrating the implementation steps of step S3. Detailed Implementation

[0016] Please see Figures 1 to 3 A control method for a robotic arm, the method comprising the following steps:

[0017] Step S1: Deploy miniature acoustic sensors at the joints of the robotic arm and monitor the acoustic waves during the joint operation to obtain joint operation monitoring acoustic waves; identify and label the joint operation monitoring acoustic waves with friction acoustic waves to output the frequency domain structure of friction acoustic waves.

[0018] Step S2: Based on the frequency domain structure simulation of the friction sound wave, derive the radial displacement index of the robotic arm joint during operation, and then perform spatial interpolation processing of the force drift of the robotic arm joint to output force drift regression data; quantify the running accuracy deviation of the joint based on the force drift regression data to obtain accuracy deviation quantification data.

[0019] Step S3: Based on the precision deviation quantification data, perform adaptive training and correction of deviations during the operation of the robotic arm to obtain adaptive correction data. Then, design the deviation correction architecture to build the deviation correction architecture and send the deviation correction architecture to the terminal to execute the control of the robotic arm.

[0020] In this embodiment of the invention, reference is made to Figure 1 The above is a flowchart illustrating the steps of a control method for a robotic arm according to the present invention. In this example, the control method for the robotic arm includes the following steps:

[0021] Step S1: Deploy miniature acoustic sensors at the joints of the robotic arm and monitor the acoustic waves during the joint operation to obtain joint operation monitoring acoustic waves; identify and label the joint operation monitoring acoustic waves with friction acoustic waves to output the frequency domain structure of friction acoustic waves.

[0022] In this embodiment of the invention, a miniature acoustic sensor is deployed at the joint of the robotic arm. The sensor has a sampling frequency of 48kHz, a sensitivity of -42dB, and a frequency response range of 20Hz to 20kHz. The sensor is fixed to the outer shell of the second joint of the robotic arm, 10mm from the bearing, using high-strength epoxy resin to ensure tight contact between the sensor and the joint surface. The sensor is connected to a data acquisition module via a shielded cable. The acquisition module is set to a sampling rate of 48kHz, a quantization accuracy of 24bit, and a acquisition duration of 120s, the time required for the robotic arm to complete one standard motion cycle. The acquired acoustic data is then processed into frames, with a frame length of 1024 sampling points and a frame shift of 512 sampling points. A Hanning window is applied to each frame of data. The function performs preprocessing to reduce spectral leakage, and then identifies and marks the framed acoustic data as frictional acoustic waves. Frictional acoustic waves are identified by setting a frequency threshold range of 1kHz to 5kHz, and a short-time energy threshold of 15dB is used for judgment. When the acoustic wave signal is within this frequency range and the energy exceeds the threshold, it is marked as a frictional acoustic wave. The marked frictional acoustic waves are subjected to a Fast Fourier Transform with 2048 transform points to obtain the spectrum of the frictional acoustic waves. Frequency components and their amplitudes are extracted from the spectrum, and a frequency-amplitude correspondence table is constructed to form the frequency domain structure of the frictional acoustic waves. This frequency domain structure contains a two-dimensional array of frequency components and their corresponding amplitudes, with a frequency resolution of 23.4Hz and an amplitude accuracy of 0.01dB. The frequency domain structure data is stored in binary format for subsequent derivation and analysis of the radial runaway index.

[0023] Step S2: Based on the frequency domain structure simulation of the friction sound wave, derive the radial displacement index of the robotic arm joint during operation, and then perform spatial interpolation processing of the force drift of the robotic arm joint to output force drift regression data; quantify the running accuracy deviation of the joint based on the force drift regression data to obtain accuracy deviation quantification data.

[0024] In this embodiment of the invention, based on the frequency domain structure of the triboelectric acoustic wave obtained in step S1, the power spectral density of the frequency components in the frequency domain structure is first calculated. The Welch method is used for estimation, with a segment length of 512 points and an overlap rate of 50%. Disorder intensity quantization is performed on the power spectral density, and the frequency narrowband sideband amplitude ratio is calculated. A center frequency of 2.5 kHz and a bandwidth of 200 Hz are selected. The amplitude ratio of the sideband to the center frequency is calculated, yielding a sideband amplitude ratio of 0.78. Subsequently, fractal feature identification of frequency amplitude abrupt changes is performed. The box-counting method is used to calculate the fractal dimension, with the box size increasing from 5 Hz to 100 Hz in a step of 5 Hz, resulting in a fractal dimension of 1.68. The peak-valley distribution density interval is calculated, with the peak detection threshold set to 1.5 times the average amplitude, yielding an average peak spacing of [value missing]. At 78Hz, with a peak height ratio variance of 0.42, the frequency disorder entropy was calculated to be 0.86. Combined with the sideband amplitude ratio, frictional acoustic wave disorder intensity data was generated, with a quantized value of 0.72. Based on this disorder intensity data and the robotic arm joint structure parameters (joint diameter 45mm, bearing inner diameter 25mm, bearing outer diameter 52mm, steel elastic modulus 210GPa, initial joint clearance 0.05mm), resonance intensity coupling simulation was performed, yielding a joint resonance coupling intensity of 0.65. Simulating the reciprocating expansion of the joint clearance, after 1000 Monte Carlo simulation iterations, the joint clearance expansion data was obtained, with an expansion amount of 0.08mm. The centrifugal vibration imbalance index was calculated to be 0.38. Integrating the center of gravity shift... With an integration time of 120 seconds and an integration step of 0.01 seconds, the integral of the center of gravity offset was obtained as 0.25 mm, and the radial slip index was derived as 0.42. Subsequently, spatial interpolation of the force drift was performed, and the time-varying characteristics of the disordered intensity data of the frictional acoustic waves were analyzed using a short-time Fourier transform with a window length of 512 points and a window shift of 256 points, resulting in a disordered time-varying intensity of 0.68. The skewness of the energy distribution density in different frequency bands was calculated: 0.56 for the low-frequency band (1kHz-2kHz), 0.72 for the mid-frequency band (2kHz-3.5kHz), and 0.48 for the high-frequency band (3.5kHz-5kHz). Based on the radial slip index, the slip frequency was calculated to be 18 Hz, the average displacement amplitude difference was 0.12 mm, and the acceleration change was calculated to be 1.54. Spatial force offset quantification was performed using m / s², with a friction contact point offset of 0.15 mm, axial acceleration components of 0.85 m / s² for the x-axis, 1.28 m / s² for the y-axis, and 0.32 m / s² for the z-axis, torsional torque offset of 0.28 N·m, and out-of-roundness simulation of 0.18 mm. This yielded spatial force offset data. Cubic spline interpolation was then used to interpolate the force drift spatially, with 200 interpolation points, resulting in force drift interpolation data. Nonlinear regression analysis was performed on the interpolation data using Gaussian process regression with a radial basis function kernel and a length scale parameter of 0.5. This output force drift regression data. Based on the regression data and the radial axial movement index, joint running accuracy deviation was quantified, with a calculated position deviation of 0.With a diameter of 22mm, an angular deviation of 0.18 degrees, and a speed deviation of 0.15mm / s, the quantitative data for accuracy deviation were obtained.

[0025] Step S3: Based on the precision deviation quantification data, perform adaptive training and correction of deviations during the operation of the robotic arm to obtain adaptive correction data. Then, design the deviation correction architecture to build the deviation correction architecture and send the deviation correction architecture to the terminal to execute the control of the robotic arm.

[0026] In this embodiment of the invention, the precision deviation quantization data obtained in step S2 is normalized using a maximum-minimum normalization method. The normalization ranges for position deviation, angle deviation, and velocity deviation are set to 0 to 1, resulting in normalized precision deviation data. The normalized values ​​for position deviation are 0.65, angle deviation, and velocity deviation are 0.58 and 0.48, respectively. Based on the normalized data, deviation increment fitting is performed over time using a polynomial fitting method with a polynomial order of 3, a fitting time window of 120 seconds, and a time step of 0.5 seconds. The coefficients of the position deviation increment fitting function are obtained as [0.0002, 0.0015, 0.0125, 0.6500]. The coefficients of the angle deviation increment fitting function are [0.0001, 0.0012, 0.0105, 0.5800], and the coefficients of the velocity deviation increment fitting function are [0.0001, 0.0010, 0.0085, 0.4800]. The precision deviation increment fitting data is output. Feature extraction is performed on the fitted data using a sliding window method with a window length of 10 seconds and a sliding step of 2 seconds. The statistical features of the deviation within each window are calculated, including the mean position deviation (0.68), variance position deviation (0.023), mean angle deviation (0.61), variance angle deviation (0.018), mean velocity deviation (0.52), and variance velocity deviation (0.015). The deviation feature data is obtained, and a multi-dimensional feature vector is constructed. The feature vector has a dimension of 12, including the mean, variance, and first and second differences of three types of biases. Principal component analysis (PCA) is used to reduce the dimensionality of the feature vector, retaining principal components with a contribution rate of 95%. After dimensionality reduction, the dimension is 5, resulting in a dimensionality-reduced feature representation. A gated recurrent unit (GRU) network is constructed based on the dimensionality-reduced features. The network structure includes 5 input layer nodes, 32 hidden layer nodes, and 3 output layer nodes. The activation function is tanh, the initial bias of the forget gate is set to 1.0, the learning rate is 0.001, the batch size is 32, the number of training epochs is 500, the training dataset size is 1000 samples, and the validation set size is 200 samples. An early stopping strategy is adopted during training; training is stopped if the loss does not decrease after 10 consecutive validation epochs, and the corrected and optimized training data is obtained. The model's root mean square error on the validation set was 0.035. Based on the trained model, adaptive training and correction of deviations during the robotic arm's operation were performed. Real-time input of current state features resulted in a position correction of -0.18 mm, an angle correction of -0.15 degrees, and a speed correction of -0.12 mm / s, yielding adaptively corrected deviation data. A random forest algorithm was used to design a deviation correction architecture for this data. Random forest parameters were set to 100 trees, a maximum depth of 15, a minimum number of leaf node samples of 5, and mean square error as the feature splitting criterion. The deviation correction architecture was constructed, comprising a feature extraction module, a feature dimensionality reduction module, a prediction module, and an execution module. This architecture was compiled into a binary file of size 2.The 5MB deviation correction architecture is sent to the terminal controller via TCP / IP protocol at a transmission rate of 10Mbps. The terminal controller is a TC-500 model with a processor clock speed of 1.2GHz and 2GB of memory. After receiving the correction architecture, the terminal loads it into memory, starts the execution thread, and uses a sampling frequency of 100Hz. The correction command latency is less than 5ms, enabling precise control of the robotic arm.

[0027] Preferably, step S1 includes the following steps:

[0028] Miniature acoustic sensors are deployed at the joints of the robotic arm to monitor the acoustic waves during the joint operation, thus obtaining the acoustic waves for joint operation monitoring.

[0029] The acoustic waves used for joint movement monitoring are processed into frames to generate monitoring frame acoustic waves.

[0030] Frictional sound waves are obtained by identifying and marking the monitored frame-by-frame sound waves.

[0031] The friction sound wave is frequency-domain converted to output the frequency domain structure of the friction sound wave.

[0032] In this embodiment of the invention, piezoelectric miniature acoustic wave sensors are fixedly installed at each rotary joint, sliding joint, and end effector connection of the robotic arm. The frequency response range of these sensors is 1kHz to 50kHz, the sensitivity is 58 dB, and the sampling frequency is set to 192000 Hz. They are connected to the synchronous acquisition module via shielded wires to ensure that the signal is not affected by electromagnetic interference during transmission. Real-time acoustic wave signals generated at the joints are acquired during the operation of the robotic arm. After entering the acquisition module, the acoustic wave signals are first preprocessed by a bandpass filter with a filtering range set to 1 kHz to 20 kHz to remove low-frequency mechanical noise and high-frequency electromagnetic interference components. The filtered signals are stored in the sampling buffer for subsequent acoustic wave framing processing. A sliding window method was used to segment the filtered continuous acoustic signal into time series segments. The window length was set to 2048 sampling points, and the window movement step size was set to 512 sampling points to form partially overlapping signal frames to ensure time continuity. Each frame of signal was weighted by a Hamming window function before fast Fourier transform to reduce edge effects. After framing, a continuous monitoring frame acoustic signal dataset was generated. Each acoustic frame contained a corresponding time index and amplitude value. The average energy, maximum peak amplitude, and frequency distribution characteristics were calculated within each time frame to generate monitoring frame acoustic signal matrix data. The data dimension was M×N, where M was the number of frames and N was the number of sampling points per frame. In an embodiment of identifying and marking friction acoustic waves in monitored framed acoustic waves, the energy envelope of each frame signal is extracted from the monitored framed acoustic wave matrix. A peak detection algorithm is used to determine the local maximum energy point. The presence of friction feature signals is determined by comparing the energy change rate of adjacent frames. When the energy change rate of adjacent frames is greater than 0.3 and exceeds 5 consecutive frames, it is determined to be a friction event. The signals within these frames are further envelope demodulated, and the carrier frequency range of the envelope signal is extracted. Friction acoustic wave features are marked in the frequency band with a frequency band energy ratio greater than 0.35, generating friction acoustic wave marking data. Each mark contains a time index, frequency band range, peak amplitude, and energy ratio parameters. This data is structured and stored as a friction acoustic wave dataset for subsequent frequency domain conversion. In the embodiment of frequency domain conversion of the friction sound wave, a fast Fourier transform is performed on each frame of the friction sound wave marker data, using an FFT of 4096 points. The converted spectrum data is then normalized to normalize the amplitude range to between 0 and 1. The main energy concentration region is extracted from the spectrum results, and the power spectral density is calculated. The frequency resolution is set to 46.9 Hz. Simultaneously, a short-time Fourier transform is used to obtain the time-frequency distribution matrix of the friction sound wave. This matrix has time as the horizontal axis, frequency as the vertical axis, and energy intensity as the matrix value, ultimately generating the frequency domain structure data of the friction sound wave. This data includes five parameters: frequency distribution, bandwidth, center frequency, energy density, and frequency change rate, which are used for subsequent radial runaway exponent derivation and force drift interpolation processing.

[0033] Preferably, step S2 includes:

[0034] The frequency domain structure of friction acoustic waves is quantized to generate disordered intensity data of friction acoustic waves.

[0035] The radial runout index during the operation of the robotic arm joint is derived by simulating the random intensity data of frictional acoustic waves.

[0036] Based on the disordered intensity data of frictional acoustic waves and the radial runout index, the force drift space interpolation of the robotic arm joint is processed to obtain the force drift interpolation data.

[0037] Nonlinear regression analysis is performed on the force drift interpolation data to output force drift regression data;

[0038] Based on the force drift regression data and radial runout index, the running accuracy deviation of the joint is quantified to obtain the accuracy deviation quantification data.

[0039] As an example of the present invention, reference is made to... Figure 2 As shown, step S2 in this example includes:

[0040] S21: Perform disordered intensity quantization on the frequency domain structure of frictional acoustic waves to generate disordered intensity data of frictional acoustic waves;

[0041] In this embodiment of the invention, in the example of disordered intensity quantization of the frequency domain structure of friction acoustic waves, the amplitude distribution curves of each frequency band in the frequency domain structure data of friction acoustic waves are first extracted. The frequency range is set to 1kHz to 20kHz, and the frequency resolution is 46.9Hz. The frequency fluctuation intensity sequence is obtained by calculating the sum of squares of the amplitude differences between adjacent frequency bands. Then, the fluctuation intensity sequence is smoothed by moving average with a sliding window length of 10 frequency points to remove high-frequency noise components. In the smoothed data, the peak-valley detection algorithm is used to extract the spacing distribution between local maximum and minimum values, and the frequency band is calculated. The amplitude ratio is used to reflect the instability characteristics of the signal. Within each frequency band, the peak height variance is calculated to measure the fluctuation complexity. Then, the Shannon entropy method is used to calculate the information entropy of the amplitude sequence. The entropy value and the amplitude ratio variance are weighted and superimposed with a weight ratio of 0.6:0.4 to generate frequency disorder entropy data. The disorder intensity index of the friction sound wave is calculated by multiplying the disorder entropy value by the frequency band energy ratio. Finally, the disorder intensity data of the friction sound wave is output. Each frequency band corresponds to a disorder intensity value, forming a disorder intensity vector of length N, where N is the number of frequency bands. This data is used as the input parameter for the subsequent radial runaway index calculation.

[0042] In another embodiment, the frequency domain structure data of the friction acoustic wave obtained in step S1 is subjected to disordered intensity quantization. First, the narrowband frequency sideband amplitude ratio of the frequency domain structure is calculated. A frequency band with a center frequency of 3.5kHz and a bandwidth of 500Hz is selected as the analysis window. The main frequency component is identified within this window, and the amplitude of the main frequency component is -28dB. Then, the amplitude of the sidebands within ±200Hz on both sides of the main frequency is calculated. The amplitude of the left sideband is -35dB, and the amplitude of the right sideband is -37dB. The amplitude ratio of the sidebands to the main frequency is calculated. The amplitude ratio of the left sideband is 0.71, and the amplitude ratio of the right sideband is 0.65. The average value of 0.68 is taken as the frequency narrowband sideband amplitude ratio. Then, the friction... The frequency-amplitude abrupt change fractal features of the acoustic wave frequency domain structure were identified. A change-point detection algorithm was used to identify abrupt changes in the spectrum, with a detection threshold of 5 dB. A total of 28 abrupt changes were detected in the 2 kHz to 10 kHz frequency band, and the distribution density of these abrupt changes was calculated to be 3.5 abrupt changes / kHz. Fractal dimension analysis was performed on these abrupt changes using box counting, with the box size increasing from 10 Hz to 500 Hz in 10 Hz increments. A double logarithmic curve was plotted, and a straight line was fitted using the least squares method. The negative value of the slope represents the fractal dimension, which was calculated to be 1.72, indicating strong fractal characteristics in the spectrum. Subsequently, the peak-valley distribution density of the amplitude abrupt change fractal features was analyzed. Peak interval calculation was performed, using a peak detection algorithm to identify peak points in the spectrum. The peak detection threshold was set to the local mean plus two standard deviations. A total of 42 peak points were detected in the 2kHz to 10kHz frequency band. The average interval of the peak points was calculated to be 190.5Hz, the standard deviation was 68.3Hz, and the coefficient of variation was 0.358. The peak height proportional variance was calculated, with an average peak height of -32.5dB, a standard deviation of 7.8dB, and a height proportional variance of 0.24. Based on the height proportional variance, the disorder entropy value of the amplitude abrupt fractal characteristics was differentiated, and the sample entropy algorithm was used to calculate the disorder of the spectrum. The embedding dimension was set to 2, and the tolerance parameter was set to 0.15 times the standard deviation. The calculated sample entropy value was 1.85, indicating that the spectrum has high complexity and disorder. The disorder intensity was quantized based on the frequency disorder entropy value and the frequency narrowband sideband amplitude ratio. A weighted fusion method was used, with an entropy weight of 0.6 and an amplitude ratio weight of 0.4. The calculated disorder intensity quantization value was 0.78, with the quantization value ranging from 0 to 1. The larger the value, the higher the degree of disorder. Finally, the disorder intensity data of the friction sound wave was generated, including the overall disorder intensity value and the disorder intensity distribution of each frequency band. The disorder intensity of the low frequency band 1kHz-3kHz was 0.65, the disorder intensity of the mid frequency band 3kHz-6kHz was 0.82, and the disorder intensity of the high frequency band 6kHz-10kHz was 0.72.

[0043] S22: The radial runout index during the operation of the robotic arm joint is derived by simulating the disordered intensity data of frictional acoustic waves.

[0044] In this embodiment of the invention, in the example of deriving the radial runout index of a robotic arm joint during operation based on disordered intensity data of frictional acoustic waves, the mass, structural stiffness, and frictional characteristics of the robotic arm joint are used as known parameters. It is assumed that the joint mass is 1.25 kg and the bearing radial stiffness is 2.3 × 10⁻⁶. 5 The N / m value is used to standardize the disordered intensity data of the frictional sound waves and treat it as the energy input signal of local frictional disturbance in the joint. The joint resonance coupling strength is derived by using the distribution of this disturbance energy input in the frequency domain. The resonance energy density of each frequency band is calculated and compared with the natural frequency range of the joint. When the resonance energy density exceeds 1.4 times the average energy density of the natural frequency range, it is marked as a resonance peak. By statistically analyzing the relationship between the resonance peak distribution and the rate of change of disordered intensity, the radial vibration amplitude increment is calculated. This increment is verified with a micron-level measuring device and set to a maximum of 0.025 mm. The radial runaway index is defined based on the average growth rate of the amplitude increment sequence. The higher the runaway index value, the stronger the radial fluctuation during joint operation. This index is used as an important physical parameter for quantifying joint runaway behavior and is input into the subsequent force drift spatial interpolation analysis.

[0045] In another embodiment, the radial runout index during the operation of the robotic arm joint is derived based on the disordered intensity data of triboelectric acoustic waves. First, the structural design parameters of the robotic arm joint are obtained, including a joint shaft diameter of 35mm, a bearing inner diameter of 35mm, a bearing outer diameter of 72mm, a bearing width of 17mm, a deep groove ball bearing type, 12 steel balls with a diameter of 11.5mm, a bearing material of GCr15 bearing steel, an elastic modulus of 210GPa, and a Poisson's ratio of 0.3, and a joint housing material of aluminum alloy with an elastic modulus of 70GPa and a Poisson's ratio of 0.33. The mass parameters of the joint are then obtained: shaft mass of 1.2kg, bearing mass of 0.35kg, and housing mass of 0.85kg. Finally, the stiffness parameters of the joint are obtained: axial stiffness of 2.8 × 10⁻⁶. 5 N / m, radial stiffness is 3.5×10 5 N / m, torsional stiffness is 4.2×10 3Based on the structural design of the robotic arm joint, the joint clearance at the factory was obtained. The shaft and bearing inner ring are interference-fitted with an interference of 0.01 mm, the bearing outer ring and housing are interference-fitted with an interference of 0.015 mm, and the bearing internal radial clearance is 0.02 mm. Resonance intensity coupling simulation of the joint's mass and stiffness was performed based on the disordered intensity data of frictional acoustic waves. A joint dynamic model was established using the finite element method, with hexahedral solid elements selected as the element type, element size 1 mm, and a total of approximately 85,000 elements. The boundary conditions were set as follows: one end of the shaft is fixed, and a torque of 15 N·m is applied to the other end at a rotational speed of 60 rpm. The natural frequencies of the joint under this condition were calculated: the first natural frequency is 125 Hz, the second natural frequency is 320 Hz, and the third natural frequency is 580 Hz. Comparison and analysis of the disordered intensity data of frictional acoustic waves with the natural frequencies revealed that the disordered intensity reaches a peak value of 0.85 near the second natural frequency, indicating the presence of resonance. The amplification factor at the resonant frequency was calculated to be 4.2, and the joint resonant coupling strength was obtained as 0.72. Based on the friction acoustic wave disorder intensity data and the joint resonant coupling strength, a random process simulation of the joint clearance reciprocating expansion at the factory was performed. The Monte Carlo method was used, with 10,000 simulations, a time step of 0.001 seconds, and a total simulation time of 300 seconds. Considering the clearance expansion caused by friction and wear, the wear rate is proportional to the disorder intensity, with a proportionality coefficient of 0.005 mm / hour / unit disorder intensity. The simulation showed that the clearance expansion after 300 hours of operation was 0.117 mm. The joint clearance expansion data was output. Based on the joint clearance expansion data and the joint resonant coupling strength, the centrifugal vibration imbalance index of the joint movement was determined. The spectral analysis method was used to calculate the ratio of vibration displacement amplitude to clearance size. This ratio was 1.85, indicating the presence of significant centrifugal vibration. The calculated peak vibration acceleration was 12.5 m / s², and the root mean square vibration velocity was 15.8 mm / s. According to ISO... According to the 10816 vibration severity standard, this vibration belongs to level C. The calculated centrifugal vibration imbalance index is 0.65, with an index range of 0 to 1. The larger the value, the higher the degree of imbalance. Based on the centrifugal vibration imbalance index, the center of gravity offset of the joint movement is integrated using a numerical integration method with an integration step size of 0.001 seconds and an integration time of 300 seconds. The calculated center of gravity offset in the x-direction is 0.085 mm, the center of gravity offset in the y-direction is 0.092 mm, the combined center of gravity offset is 0.125 mm, and the center of gravity offset integral is 37.5 mm·s. Based on the centrifugal vibration imbalance index and the center of gravity offset integral, the radial movement index during the operation of the robotic arm joint is derived. Using a weighted fusion method, the imbalance index weight is 0.7, and the offset integral weight is 0.3. The calculated radial movement index is 0.58, with an index range of 0 to 1. The larger the value, the higher the degree of movement.

[0046] S23: Based on the disordered intensity data of frictional acoustic waves and the radial displacement index, perform spatial interpolation processing of the force drift of the robotic arm joint to obtain force drift interpolation data;

[0047] In this embodiment of the invention, the embodiment of spatial interpolation processing of force drift of robotic arm joint based on disordered intensity data of frictional acoustic waves and radial runout index firstly performs joint normalization processing on the radial runout index of the joint and the disordered intensity data, keeping their numerical range between 0 and 1. Then, a joint force space model is established in a three-dimensional coordinate system, where the x, y, and z axes represent the torque distribution of the joint in three spatial directions, respectively. Using the energy change trend corresponding to the radial runout index, a discrete set of force offset points is generated in three-dimensional space. The number of points in this set is set to 500. The point set is spatially continuous through a three-dimensional spline interpolation algorithm. A smooth surface is constructed using a cubic spline function to represent the spatial distribution of force drift. In the interpolation calculation, the node spacing is controlled at 0.5 mm to ensure the smoothness of the interpolation. The force offset of each node is calculated, with the offset unit being N·m. Finally, a force drift interpolation data matrix is ​​generated. Each element in this matrix contains force components in three directions and their corresponding displacement offsets, which are used for subsequent nonlinear regression analysis.

[0048] In another embodiment, spatial interpolation processing of the force drift of the robotic arm joint is performed based on the disordered intensity data of friction acoustic waves and the radial runout index. First, the time-varying characteristics of the disordered intensity data of friction acoustic waves are analyzed using the short-time Fourier transform method with a window length of 1024 points and a window shift of 512 points. The time-varying spectrum is calculated, and the characteristics of the spectrum changing with time are analyzed. The disordered intensity value at each time point is extracted to form a disordered intensity time series. The time average value of the disordered intensity is calculated to be 0.78, the standard deviation is 0.12, and the coefficient of variation is 0.154. The disordered time-varying intensity data is output. Based on the disordered time-varying intensity, the energy distribution density skewness of the disordered intensity data of friction acoustic waves in different frequency bands is evaluated, and the spectrum is divided into 5 frequency bands: 1kHz-2kHz. The energy distribution was calculated in the frequency bands 2kHz-4kHz, 4kHz-6kHz, 6kHz-8kHz, and 8kHz-10kHz, accounting for 12%, 28%, 35%, 18%, and 7% of the total energy, respectively. The skewness of the energy distribution was calculated by dividing the third central moment by the cube of the standard deviation, yielding a skewness value of 0.65, indicating a bias towards higher frequencies. The skewness of the energy distribution density in each frequency band was calculated to be 0.32, 0.58, 0.75, 0.48, and 0.25, respectively. Based on the radial runaway index, the runaway frequency of the robotic arm joints was counted using the zero-crossing rate method, with a threshold set to 1.2 times the average value. The calculated runaway frequency was 22.5Hz. The traversal cycle is 0.044 seconds, and approximately 4100 traversals occur within one standard work cycle. The radial displacement amplitude is analyzed, and the peak displacement is extracted using the peak detection method. The average peak value is 0.085 mm, the maximum value is 0.132 mm, the minimum value is 0.042 mm, and the standard deviation is 0.023 mm. The mean difference of displacement amplitude is calculated to be 0.09 mm. Based on the mean difference of displacement amplitude, the acceleration change during the traversal process of the robotic arm joint is calculated. The acceleration is calculated using the central difference method. The average acceleration is 8.5 m / s², the maximum value is 15.2 m / s², the minimum value is 3.8 m / s², and the standard deviation is 2.8 m / s². The spatial distribution of the robotic arm joint is analyzed based on the energy distribution density skewness data and the acceleration change. Force offset quantification: First, the offset of the friction contact point is analyzed based on the energy distribution density skewness data, establishing a mapping relationship between skewness and contact point offset with a mapping coefficient of 0.15 mm / unit skewness. The calculated offset of the friction contact point is 0.098 mm. The acceleration variation is then analyzed by analyzing the directional component characteristics, decomposing it into radial and tangential components. The radial acceleration component is 7.2 m / s², and the tangential acceleration component is 4.5 m / s². The calculated acceleration direction angle is 32°, obtaining the axial acceleration direction component data. Based on the friction contact point offset and the axial acceleration direction component, torsional moment offset fitting is performed using a polynomial fitting method with a fitting order of 3 and fitting coefficients of [0.0025, 0.018, 0.12, 0.018].

[35] The torsional torque offset was calculated to be 0.42 N·m, and the torsional torque offset rate was 2.8%. Based on the torsional torque offset data and the directional component of the axial acceleration, a numerical simulation of the out-of-roundness of the minimum circumscribed circle of the circular motion was performed. The least squares method was used to fit the circle, and the maximum deviation from the actual trajectory point to the fitted circle was calculated to be 0.128 mm, the average deviation to be 0.075 mm, and the standard deviation to be 0.032 mm. The out-of-roundness was calculated to be 0.128 mm, and the simulated out-of-roundness value was obtained. Based on the torsional torque offset data and the simulated out-of-roundness value, a mechanical... Spatial force offset quantification of the robotic arm joint was performed, establishing a mapping relationship between torque offset and spatial force with a mapping coefficient of 2.5 N / N·m. The calculated spatial force offset was 1.05 N, with the direction of the spatial force offset aligned with the acceleration direction. This spatial force offset data was then processed using spatial force drift interpolation. Radial basis function interpolation was employed, with multiple quadratic functions selected as the basis functions. The smoothing parameter was set to 0.01, and 200 interpolation points were used, covering the entire workspace of the joint motion. The interpolation accuracy was 0.01 N, yielding the force drift interpolation data.

[0049] S24: Perform nonlinear regression analysis on the force drift interpolation data to output force drift regression data;

[0050] In this embodiment of the invention, in the example of performing nonlinear regression analysis on force drift interpolation data, regression analysis is performed on the force drift interpolation data matrix. The input variable is the force components in three-dimensional space, and the output variable is the drift displacement in each direction. A quadratic polynomial nonlinear regression algorithm is used for fitting. The coefficients of each term in the regression equation are determined by the least squares method. During the data fitting process, the sum of squared residuals is minimized to improve the fitting accuracy. The regression coefficients are repeatedly calculated for each joint position point and cross-validation is performed. The number of validation samples is set to 100 groups. After the validation is completed, the force drift regression data is output. This data is structured into a vector form, and each component corresponds to the predicted displacement regression value of a force point, with the unit being mm.

[0051] In another embodiment, nonlinear regression analysis is performed on the obtained force drift interpolation data. First, the interpolation data is preprocessed using a moving average filter to remove noise. The window length is 5 points, which improves the smoothness of the filtered data by 25%. Outlier detection is then performed on the filtered data using the 3σ criterion, identifying 12 outlier points, accounting for 2.4% of the total points. These outliers are replaced with local means to ensure data continuity and smoothness. Feature extraction is then performed on the processed data. Time-domain features include mean, standard deviation, gamma, skewness, and kurtosis. Frequency-domain features include dominant frequency components and frequency band energy distribution. A feature vector with a dimension of 15 is constructed. The extracted feature vectors were used to construct a nonlinear regression model. A Gaussian process regression method was employed, with a radial basis function as the kernel function. The kernel parameter length scale was set to 0.5, the signal variance to 1.0, and the noise variance to 0.01. The training set size was 400 points, and the test set size was 100 points. Five-fold cross-validation was used to evaluate the model performance. The root mean square error was 0.025 mm, and the coefficient of determination (R²) was 0.92, indicating a good model fit. The trained Gaussian process regression model was used to predict force drift interpolation data, yielding force drift regression data. The regression data dimension was 500×3, representing the three-dimensional force drift prediction values ​​for 500 spatial points. The x-axis component of the force drift regression data... The mean was 0.15 mm, and the standard deviation was 0.04 mm. The mean of the y-axis component was 0.18 mm, and the standard deviation was 0.05 mm. The mean of the z-axis component was 0.08 mm, and the standard deviation was 0.02 mm. Uncertainty quantification was performed on the regression data, and the 95% confidence interval for each prediction point was calculated. The average confidence interval width was 0.06 mm, indicating high reliability of the prediction results. A force drift regression data structure was constructed, including regression prediction values, prediction uncertainty, model parameters, etc., providing basic data for subsequent quantification of joint operation accuracy deviation. Furthermore, sensitivity analysis was performed on the regression model, observing output changes through perturbation input characteristics to identify key influencing factors and discover x. The force drift in the x-axis and y-axis directions has the greatest impact on the model output, with sensitivity coefficients of 0.65 and 0.72, respectively. The influence in the z-axis direction is relatively small, with a sensitivity coefficient of 0.28. Based on the sensitivity analysis results, the regression model is optimized by adjusting the feature weights. The feature weights of the x-axis and y-axis are increased by 1.2 times, while the feature weight of the z-axis is decreased by 0.8 times. After optimization, the root mean square error of the model is reduced to 0.022 mm, and the coefficient of determination R² is increased to 0.94. The final output is the optimized force drift regression data, which includes timestamps, spatial coordinates, three-axis force drift values ​​and their uncertainty estimates. The data sampling rate is 100 Hz, the time span is 180 seconds, and there are a total of 18,000 data points.

[0052] S25: Based on the force drift regression data and radial displacement index, the running accuracy deviation of the joint is quantified to obtain the accuracy deviation quantification data.

[0053] In this embodiment of the invention, in the embodiment of quantifying the joint's operational accuracy deviation based on the force drift regression data and radial runout index, the force drift regression data and the radial runout index are jointly analyzed. Synchronous interpolation is performed on both on the time axis to ensure time correspondence. The difference between the radial runout index change rate and the force drift regression value at the same time point is normalized. The deviation value is obtained by calculating the Euclidean distance between the regression predicted displacement and the actual joint angle sensor measured displacement. The root mean square operation is performed on the deviation values ​​over consecutive time periods to obtain the overall operational accuracy deviation index of the joint. The deviation results are output as accuracy deviation quantification data in the form of a time series, with the unit being mm. This data reflects the displacement error distribution of the joint throughout the entire motion cycle and is used in the subsequent adaptive correction stage of the deviation.

[0054] In another embodiment, based on the output force drift regression data and the derived radial runout index, a kinematic model of the robotic arm joints is first established. Using the Denavit-Hartenberg parameter method, the DH parameters for joint 1 are: a1=0mm, d1=150mm, α1=90°, θ1 is the joint variable; the DH parameters for joint 2 are: a2=250mm, d2=0mm, α2=0°, θ2 is the joint variable; and the DH parameters for joint 3 are: a3=220mm, d3=0mm, α3=0°, θ3 is the joint variable. Based on the DH parameters, a joint coordinate system transformation matrix is ​​constructed, and the mapping relationship between the position of the robotic arm end effector and the joint angles is calculated. A Jacobi matrix was established to represent the relationship between the joint angle error and the end-effector position error. The Jacobi matrix has a dimension of 6×6 and represents the relationship between the end-effector errors of the six degrees of freedom and the six joint angle errors. Force drift regression data were converted into joint angle deviations. Using inverse kinematics, the joint angle deviation corresponding to each force drift point was calculated. The mean angle deviation for joint 1 was 0.12°, and the standard deviation was 0.04°; the mean angle deviation for joint 2 was 0.18°, and the standard deviation was 0.05°; the mean angle deviation for joint 3 was 0.15°, and the standard deviation was 0.04°. The joint angle deviations were corrected using a radial runout index, with a correction coefficient of 1 + 0.5 × radial runout index, where the radial runout index was 0. With a correction factor of 1.325, the mean angle deviation of joint 1 is 0.16°, the mean angle deviation of joint 2 is 0.24°, and the mean angle deviation of joint 3 is 0.20°. Based on the corrected joint angle deviations, the end effector position deviation is calculated using the forward kinematics method. The calculated end effector position deviations are as follows: x-axis component: mean 0.85 mm, standard deviation 0.22 mm; y-axis component: mean 0.92 mm, standard deviation 0.25 mm; z-axis component: mean 0.65 mm, standard deviation 0.18 mm; roll component: mean 0.28°, standard deviation 0.08°; pitch component: mean 0.32°, standard deviation 0.08°. The deviation was 0.09°, the mean of the yaw component was 0.25°, and the standard deviation was 0.07°. A comprehensive evaluation of position and attitude deviations was conducted using the weighted Euclidean distance method, with a position weight of 0.7 and an attitude weight of 0.3. The calculated comprehensive deviation index was 0.82. Normalizing the comprehensive deviation index to the range of 0 to 1 yielded a normalized comprehensive deviation of 0.68. Based on the normalized comprehensive deviation, a precision deviation quantification model was constructed. A piecewise linear mapping method was used to map the normalized comprehensive deviation to a precision level, which was divided into five levels: Excellent (0-0.2), Good (0.2-0.4), Average (0.4-0.6), Poor (0.6-0.8), and Inferior (0.8-1).[0] The current accuracy level is "poor". Time series analysis is performed on the accuracy deviation using an autoregressive integral moving average (ARIMA) model with order parameters p=2, d=1, and q=1 to predict the trend of accuracy deviation over the next 30 seconds. The prediction results show that the accuracy deviation is increasing and is expected to reach 0.75 after 30 seconds. A quantitative data structure for accuracy deviation is constructed, including information such as joint angle deviation, end-effector position deviation, attitude deviation, comprehensive deviation index, accuracy level, and time series prediction results. The data sampling rate is 100Hz, the time span is 180 seconds, and there are a total of 18,000 data points. Each data point includes a timestamp, three joint angle deviations, six-dimensional end-effector deviation, and a comprehensive deviation index. In addition to labeling and accuracy levels, spatial distribution analysis was performed on the quantified accuracy deviation data. An accuracy deviation distribution map was drawn within the workspace, which was divided into a 10×10×10 grid. Local accuracy deviation was calculated for each grid point. It was found that the accuracy deviation was larger at the edge of the robotic arm's workspace, reaching a maximum of 0.95, while the accuracy deviation was smaller in the central area, with a minimum of 0.42. Based on the spatial distribution characteristics, the quantified accuracy deviation data was partitioned into three regions: high-precision, medium-precision, and low-precision. This partitioning strategy provides a basis for subsequent adaptive training and correction of the deviation. Finally, the quantified accuracy deviation data was output, providing a data foundation for the adaptive training and correction of the deviation in step S3.

[0055] Preferably, performing disordered intensity quantization includes:

[0056] The narrowband frequency sideband amplitude ratio of the triboacoustic wave frequency domain structure is calculated to obtain the frequency narrowband sideband amplitude ratio.

[0057] Based on the frequency narrowband sideband amplitude ratio, frequency amplitude abrupt fractal features of the frequency domain structure of the friction acoustic wave are identified, and amplitude abrupt fractal features are obtained.

[0058] The peak-valley distribution density interval is calculated for the fractal characteristics of amplitude abrupt change, and then the height ratio variance of the peaks is calculated.

[0059] Based on the high proportional variance, the disordered entropy value of the amplitude abrupt fractal feature is differentiated by the disordered entropy value to obtain the frequency disordered entropy value.

[0060] Disorder intensity is quantized based on the frequency disorder entropy value and the frequency narrowband sideband amplitude ratio to generate triboacoustic wave disorder intensity data.

[0061] In this embodiment of the invention, the narrowband frequency sideband amplitude ratio is calculated for the triboelectric acoustic wave frequency domain structure obtained in step S1. First, the main frequency components are extracted from the frequency domain structure. A peak detection algorithm is used, with the detection threshold set to 1.5 times the average amplitude and the minimum peak spacing set to 200Hz. A total of 15 main frequency peaks are detected, mainly distributed in the range of 2kHz to 10kHz. The five frequency points with the highest amplitude are selected as the center frequencies, namely 2.8kHz, 4.2kHz, 5.6kHz, 7.3kHz, and 8.9kHz. For each center frequency, a narrowband analysis bandwidth of 200Hz is set, which is the spectrum within a 100Hz range to the left and right of the center frequency. Spectral data within each narrowband is extracted, including frequency points and their corresponding amplitudes. For a narrowband with a center frequency of 2.8kHz, the frequency range is 2.7kHz to 2.9kHz, with 9 sampling points and a frequency resolution of 23.4Hz. The average amplitude of the left and right sidebands within the narrowband is calculated. The average amplitude of the left sideband (2.7kHz-2.8kHz) is -32.5dB, and the average amplitude of the right sideband (2.8kHz-2.9kHz) is... The amplitude is -33.8dB, and the amplitude at the center frequency is -28.2dB. The ratio of the sideband amplitude to the center frequency amplitude is calculated. The left sideband amplitude ratio is 0.87, and the right sideband amplitude ratio is 0.83. The average of the left and right sideband amplitude ratios is taken as the narrowband sideband amplitude ratio at this center frequency, resulting in a narrowband sideband amplitude ratio of 0.85 at 2.8kHz. The same method is used to calculate the narrowband sideband amplitude ratios for the other four center frequencies: 0.78 at 4.2kHz, 0.72 at 5.6kHz, and 0.75 at 7.3kHz. The sideband amplitude ratio is 0.68, and the narrowband amplitude ratio at 8.9kHz is 0.65. The weighted average of the narrowband amplitude ratios at the five center frequencies is calculated, with the weights proportional to the amplitude of each center frequency, resulting in a frequency narrowband amplitude ratio of 0.75. This value reflects the sharpness of the triboelectric sound wave spectrum. The smaller the sideband amplitude ratio, the sharper the spectral peaks and the more irregular the spectral structure. A frequency narrowband amplitude ratio data structure is constructed, containing five center frequency points, their respective left and right sideband amplitude ratios, and the weighted average, providing basic data for subsequent identification of frequency amplitude abrupt fractal features.

[0062] Based on the calculated frequency narrowband sideband amplitude ratio, the frequency amplitude abrupt change fractal feature identification of the frequency domain structure of the friction acoustic wave was performed. First, the frequency domain structure data was preprocessed, using median filtering to remove isolated noise points with a filter window length of 5 points. The filtered spectrum data was then normalized, adjusting the amplitude range to between 0 and 1 for easier subsequent analysis. The CUSUM (Cumulative Sum Control Chart) algorithm was used to detect abrupt changes in the spectrum. The algorithm parameters were set as follows: detection threshold h = 0.05, drift parameter k = 0.01, and confidence level of 95%. A total of 32 abrupt changes were detected in the frequency range of 2kHz to 10kHz. The mutation point frequencies were categorized as follows: 2.35kHz, 2.82kHz, 3.15kHz, 3.68kHz, 4.12kHz, 4.35kHz, 4.78kHz, 5.23kHz, 5.65kHz, 5.92kHz, 6.28kHz, 6.55kHz, 6.87kHz, 7.24kHz, 7.58kHz, 7.85kHz, 8.12kHz, 8.45kHz, 8.78kHz, 9.15kHz, 9.42kHz, and 9.85kHz. Cluster analysis was performed on the detected mutation points using a hierarchical clustering algorithm, with distance metric selected as... Using Euclidean distance and the Ward method for clustering, with 5 clusters, 5 mutation point clusters were obtained, located in the frequency ranges of 2.3-3.2kHz, 3.6-4.4kHz, 5.2-6.0kHz, 6.8-7.6kHz, and 8.1-9.5kHz, respectively. The amplitude change of each mutation point, i.e., the difference in amplitude before and after the mutation, was calculated. The amplitude change ranged from 5.2dB to 18.6dB, with an average change of 12.3dB. Fractal analysis was performed on the spatial distribution of the mutation points, and the fractal dimension was calculated using box counting, with the box size increasing from 10Hz to 500Hz in 10Hz increments. The calculated fractal dimension is 1.68, indicating that the distribution of mutation points has a certain degree of self-similarity. An amplitude mutation feature descriptor is constructed by combining the frequency narrowband sideband amplitude ratio. The feature descriptor includes: 32 mutation points, a mutation point fractal dimension of 1.68, an average mutation amplitude of 12.3 dB, 5 mutation point clusters, and a narrowband sideband amplitude ratio of 0.75. Based on these feature parameters, an amplitude mutation fractal feature vector is constructed with a dimension of 15, which includes the above 5 global features and 10 local features (the number of mutation points and the average mutation amplitude of each cluster). The amplitude mutation fractal feature data is obtained, which provides basic data for subsequent calculation of peak and valley distribution density intervals.

[0063] The peak-valley distribution density interval was calculated for the fractal characteristics of amplitude mutations. First, peak detection was performed on the frequency domain structure of the frictional acoustic wave using a local maximum detection algorithm. Detection parameters were set as follows: minimum peak height 1.5 times the average amplitude, minimum peak spacing 100Hz. A total of 45 peak points were detected within the frequency range of 2kHz to 10kHz, with peak frequencies distributed between 2.2kHz and 9.8kHz. The frequency interval between adjacent peak points was calculated, yielding 44 interval values ​​ranging from 105Hz to 485Hz, with an average interval of 175Hz and a standard deviation of 85Hz. Valley pairing was then performed on the peak points, with the lowest points on either side of each peak point serving as the corresponding valley points. Ninety valley points were detected. The peak-valley height difference was calculated, which is the difference between the amplitude of the peak point and the amplitude of the corresponding valley point. The height difference ranged from 6.5 dB to 22.8 dB, with an average height difference of 14.2 dB and a standard deviation of 4.8 dB. The peak-valley distribution density was calculated, defined as the number of peak-valley pairs per unit frequency range. The frequency range was 2 kHz to 10 kHz, with a total bandwidth of 8 kHz. There were 45 peak-valley pairs, resulting in a peak-valley distribution density of 5.63 pairs / kHz. The peak-valley distribution was segmented, dividing the 2 kHz to 10 kHz frequency range into eight equal sub-bands, each with a width of 1 kHz. The number and distribution density of peak-valley pairs within each sub-band were calculated. Subband 1 (2-3kHz): 7 peak-valley pairs, distribution density 7 pairs / kHz; Subband 2 (3-4kHz): 8 peak-valley pairs, distribution density 8 pairs / kHz; Subband 3 (4-5kHz): 9 peak-valley pairs, distribution density 9 pairs / kHz; Subband 4 (5-6kHz): 6 peak-valley pairs, distribution density 6 pairs / kHz; Subband 5 (6-7kHz): 5 peak-valley pairs, distribution density 5 pairs / kHz; Subband 6 (7-8kHz): 4 peak-valley pairs, distribution density 4 pairs / kHz; Subband 7 (8-9kHz): 3 peak-valley pairs, distribution density 3 pairs / kHz. Subband 8 (9-10kHz): There are 3 peak-valley pairs, with a distribution density of 3 pairs / kHz. The height ratio variance of the peaks is calculated. First, the height ratio of each peak is calculated, defined as the ratio of the peak height to the average peak height in the subband. 45 height ratio values ​​are obtained, ranging from 0.65 to 1.85, with an average height ratio of 1.0. The variance of the height ratio is calculated, and the variance is 0.12. This value reflects the non-uniformity of the peak height distribution. The larger the variance, the more non-uniform the peak height distribution and the more irregular the spectral structure. A peak-valley distribution feature data structure is constructed, including information such as the number of peak-valley pairs, average interval, peak-valley distribution density, and height ratio variance.

[0064] Based on the calculated height proportional variance, disordered entropy differentiation is performed on the fractal characteristics of amplitude mutations. First, an amplitude mutation sequence is constructed, arranging the mutation points in the frequency domain structure in frequency order to form a mutation amplitude sequence with a length of 32. The mutation amplitude sequence is normalized, adjusting the amplitude range to between 0 and 1 to facilitate subsequent entropy calculation. The sample entropy algorithm is used to calculate the sequence complexity, with the algorithm parameters set as follows: embedding dimension m=2, tolerance r=0.15×standard deviation. The calculated sample entropy value is 0.85, indicating that the sequence has high complexity and irregularity. Entropy analysis is performed on the frequency distribution of mutation points using the information entropy method. The frequency range of 2kHz to 10kHz is divided into 20 bins, and the probability distribution of mutation points within each bin is calculated, yielding a frequency distribution entropy value of 2.65. The same method is used to analyze the entropy distribution of mutation points. The amplitude range was divided into 20 bins, and the amplitude distribution entropy was calculated to be 2.42. The disordered entropy was differentiated using the height proportion variance, and a weighted differentiation method was employed to calculate the first derivative of the disordered entropy. With a derivative step size of 0.01, the derivative of the disordered entropy was found to be 0.28, indicating the sensitivity of the disordered entropy to the height proportion variance; a larger derivative value indicates higher sensitivity. A comprehensive disordered entropy was constructed based on sample entropy, frequency distribution entropy, amplitude distribution entropy, and their derivatives. A weighted average method was used, with the weights allocated as follows: sample entropy 0.4, frequency distribution entropy 0.3, amplitude distribution entropy 0.2, and entropy derivative 0.1. The calculated frequency disordered entropy was 1.68. The disordered entropy was normalized to the range of 0 to 1, resulting in a normalized frequency disordered entropy of 0.84. A data structure for the frequency disordered entropy was constructed, containing information such as sample entropy, frequency distribution entropy, amplitude distribution entropy, entropy derivative, and comprehensive disordered entropy.

[0065] Based on the obtained frequency disorder entropy value and the calculated frequency narrowband sideband amplitude ratio, disorder intensity quantification is performed. First, a disorder intensity assessment index system is constructed, including two primary indicators: frequency disorder entropy value and frequency narrowband sideband amplitude ratio; four secondary indicators: sample entropy, frequency distribution entropy, amplitude distribution entropy, and sideband amplitude ratio; and eight tertiary indicators: the values ​​and rates of change of each secondary indicator. A hierarchical analysis (AHP) is performed on the index system to determine the weights of each indicator. The weights for the primary indicators are: frequency disorder entropy value 0.65, frequency narrowband sideband amplitude ratio 0.35; the weights for the secondary indicators are: sample entropy 0.3, frequency distribution entropy 0.2, amplitude distribution entropy 0.15, and sideband amplitude ratio 0.35; and the weights for the tertiary indicators are equally distributed according to the weights of the secondary indicators. Based on the determined index weights, a disorder intensity quantification model is constructed, and a weighted summation method is used to calculate... The disorder intensity value is calculated to be 0.81, with a frequency disorder entropy value of 0.84 and a frequency narrowband sideband amplitude ratio of 0.75. A nonlinear mapping is then performed on this initial disorder intensity value using a sigmoid function with parameters set to a midpoint of 0.5 and a slope factor of 10, mapping the disorder intensity value to the range of 0 to 1, resulting in a mapped disorder intensity value of 0.88. Time-series smoothing is then applied to the disorder intensity value using exponential smoothing with a smoothing coefficient α = 0.3. The disorder intensity values ​​for the first five frames are 0.88, 0.85, 0.82, 0.86, and 0.84, respectively. After smoothing, the disorder intensity value for the current frame is 0.85. A data structure for the disorder intensity of the friction acoustic wave is constructed, containing information such as the disorder intensity value, frequency disorder entropy value, frequency narrowband sideband amplitude ratio, values ​​of various indices, and their weights.

[0066] Preferably, deriving the radial runout index during the operation of the robotic arm joint includes:

[0067] Obtain the structural design, mass, and stiffness of the robotic arm joints; obtain the joint clearance at the time of manufacture based on the structural design of the robotic arm joints;

[0068] Based on the disordered intensity data of frictional acoustic waves, the mass and stiffness of the joint are simulated by resonance intensity coupling to obtain the joint resonance coupling intensity.

[0069] Based on the aforementioned friction acoustic wave disorder intensity data and joint resonance coupling intensity, a random process simulation of joint gap reciprocating expansion at the factory is performed, and joint gap expansion data is output.

[0070] The eccentric vibration imbalance index of joint motion is determined based on joint space enlargement data and joint resonance coupling strength.

[0071] The center of gravity offset integral is obtained by integrating the center of gravity offset during joint movement based on the centrifugal vibration imbalance index; the radial movement index during the operation of the robotic arm joint is derived based on the centrifugal vibration imbalance index and the center of gravity offset integral.

[0072] In this embodiment of the invention, the structural design parameters of the robotic arm joints are obtained. In this embodiment, the robotic arm is a 6-axis industrial robotic arm. Joint 1 uses a harmonic reducer of model SHG-32-100-2UH with a reduction ratio of 100:1 and a rated torque of 137 N·m. The joint bearing uses an angular contact ball bearing of model 7208C with an inner diameter of 40 mm and an outer diameter of 80 mm. The joint housing material is aluminum alloy 6061-T6 with an elastic modulus of 68.9 GPa. The mass of each joint is measured: the total mass of joint 1 is 4.8 kg, the total mass of joint 2 is 3.2 kg, and the total mass of joint 3 is 2.5 kg. The stiffness characteristics of each joint are measured: the torsional stiffness of joint 1 is 2.8 × 10⁻⁶. 5 N·m / rad, the torsional stiffness of joint 2 is 1.5×10 N·m / rad. 5 N·m / rad, the torsional stiffness of joint 3 is 9.2×10 N·m / rad. 4 The natural frequencies of each joint were measured using the hammer impact method (N·m / rad). The first natural frequency of joint 1 was 285Hz. Based on the structural design of the robotic arm joints, the joint clearances at the factory were obtained. The backlash of the SHG-32-100-2UH reducer at the factory was 0.5 arcminutes, and the radial clearance of the 7208C bearing at the factory was 0.015mm. Combining the backlash of the reducer and the bearing clearance, the overall clearance of joint 1 at the factory was calculated to be 0.025mm.

[0073] Based on the disordered intensity data of the frictional acoustic waves obtained in step S4, a resonant intensity coupling simulation of the joint's mass and stiffness is performed. First, a joint dynamic model is established. Using the lumped parameter method, the joint is simplified into a mass-spring-damping system. The equivalent mass of joint 1 is 4.8 kg, and the equivalent torsional stiffness is 2.8 × 10⁻⁶. 5 With an equivalent damping ratio of 0.05 and a natural frequency of 285 Hz for joint 1, the friction acoustic wave disorder intensity data was converted into an excitation force spectrum. Using deconvolution, the disorder intensity value of 0.85 was mapped to the excitation force spectrum, with a frequency range of 0 Hz to 10 kHz. The main energy is concentrated in the range of 2 kHz to 8 kHz, and the peak value of the excitation force spectrum is 0.5 N. The excitation force spectrum was applied to the joint dynamics model, and the frequency response function (FRF) of the joint was calculated using frequency domain analysis. The FRF of joint 1 is at 28 Hz. A distinct resonance peak with an amplitude of 15 dB is observed at 5 Hz. The convolution of the excitation force spectrum and the joint FRF is calculated to obtain the joint response spectrum. The response spectrum of joint 1 has a peak value of 7.5 N at 285 Hz. The correlation between the response spectrum and the friction acoustic wave spectrum is analyzed. The cross-correlation analysis method is used to calculate the cross-correlation coefficient between the two. The cross-correlation coefficient of joint 1 is 0.72. The resonance coupling strength is calculated based on the cross-correlation coefficient. Normalization is used to map the cross-correlation coefficient to the range of 0 to 1. The resonance coupling strength of joint 1 is 0.72.

[0074] Based on the obtained triboacoustic wave disorder intensity data and the calculated joint resonance coupling strength, a stochastic process simulation of the reciprocating expansion of the joint clearance at the factory setting is performed. First, a physical model of the joint clearance expansion is established, considering three main factors: frictional wear, thermal expansion, and mechanical impact. The clearance expansion caused by frictional wear is proportional to the triboacoustic wave disorder intensity, with a proportionality coefficient of 0.05 mm. A stochastic process model of the joint clearance expansion is then established using the Markov chain Monte Carlo method. The state space represents the joint clearance value, with the initial state being the factory clearance value (0.025 mm for joint 1). The state transition probability matrix is ​​constructed based on the triboacoustic wave disorder intensity data and the resonance coupling strength. The transition probability is related to the disorder intensity... The value is proportional to the resonant coupling strength. The disordered strength value is 0.85, the resonant coupling strength of joint 1 is 0.72, and the state transition probability is 0.612. The simulation steps are set to 10,000 steps, each step corresponds to an actual time of 0.1 seconds, and the total simulation time is 1,000 seconds. In each simulation step, the joint gap is determined according to the state transition probability. If it is expanded, the expansion amount follows a log-normal distribution with distribution parameters μ=-5 and σ=0.8. The corresponding average expansion amount is 0.0068mm. After 10,000 simulation steps, the change in joint gap is statistically analyzed. The gap of joint 1 expands from the initial value of 0.025mm to the final value of 0.082mm, an increase of 0.057mm.

[0075] Based on the output joint clearance expansion data and the obtained joint resonance coupling strength, the eccentric vibration imbalance index of the joint motion is determined. First, a joint eccentric vibration model is established. Using rotor dynamics theory, the joint is simplified into a rigid rotor-elastic support system. The rotor mass of joint 1 is 4.8 kg, and the support stiffness is 2.8 × 10⁻⁶. 5The joint clearance is 0.082 mm after expansion. The unbalance of the joint is calculated using the mass eccentricity method. The unbalance is equal to the product of the mass and the eccentricity. The eccentricity is proportional to the joint clearance, with a proportionality coefficient of 0.8. The eccentricity of joint 1 is 0.066 mm, and the unbalance is 0.317 kg·mm. The centrifugal force of the joint at different speeds is calculated. The centrifugal force is equal to the product of the unbalance and the square of the angular velocity. The centrifugal force of joint 1 at the rated speed of 120 rpm (2 Hz) is 0.132 N, and at the maximum speed of 240 rpm... The centrifugal force at 4 Hz is 0.528 N. The relationship between centrifugal force and joint resonance coupling strength is analyzed using linear regression to establish a mapping relationship. The regression coefficient is 0.85, and the coefficient of determination R² is 0.92, indicating a strong correlation between the two. The resonance coupling strength of joint 1 is 0.72, and the corresponding centrifugal force correction coefficient is 1.612. The corrected centrifugal force is 0.213 N at rated speed and 0.851 N at maximum speed. The critical speed of the joint is calculated; the critical speed is equal to the natural frequency divided by 6. 0. The first critical speed of joint 1 is 285Hz / 60 = 4.75Hz, corresponding to 285rpm. Analyzing the relationship between the joint's operating speed and critical speed, the speed ratio is defined as the ratio of the operating speed to the critical speed. The speed ratio of joint 1 at rated speed is 120 / 285 = 0.421, and at maximum speed is 240 / 285 = 0.842. Based on the speed ratio and the corrected centrifugal force, the centrifugal vibration imbalance index is calculated using the weighted product method. The centrifugal vibration imbalance index is equal to the product of the speed ratio and the corrected centrifugal force multiplied by the weighting coefficient 0. 0.5, the centrifugal vibration imbalance index of joint 1 at rated speed is 0.421×0.213×0.5=0.045, and the centrifugal vibration imbalance index at maximum speed is 0.842×0.851×0.5=0.358. Calculate the average centrifugal vibration imbalance index of the joint in the actual working cycle using the time-weighted average method. The working time of joint 1 at rated speed accounts for 70%, and the working time at maximum speed accounts for 30%. The average centrifugal vibration imbalance index is 0.045×0.7+0.358×0.3=0.139.

[0076] Based on the calculated centrifugal vibration imbalance index, the integral of the center of gravity shift during joint motion is performed. First, a joint center of gravity shift model is established, using the elastic deformation theory. Under the action of centrifugal force, the joint undergoes elastic deformation, leading to a shift in the center of gravity. The amount of center of gravity shift is directly proportional to the centrifugal force and inversely proportional to the joint stiffness. The stiffness of joint 1 is 2.8 × 10⁻⁶. 5Given a centrifugal vibration imbalance index of 0.139 (N·m / rad), the instantaneous center of gravity offset of the joint is calculated using a linear mapping method, mapping the centrifugal vibration imbalance index to the center of gravity offset with a mapping coefficient of 0.5 mm. The instantaneous center of gravity offset of joint 1 is 0.139 × 0.5 = 0.0695 mm. Analyzing the directional characteristics of the center of gravity offset, a phase analysis method is used to decompose the offset into radial and tangential components. The radial component is in the same direction as the rotation radius, and the tangential component is in the same direction as the rotation radius. With the radius perpendicular, the radial component of the center of gravity offset of joint 1 accounts for 85%, and the tangential component accounts for 15%. The radial component is 0.0695 × 0.85 = 0.059 mm, and the tangential component is 0.0695 × 0.15 = 0.01 mm. The time-varying characteristics of the center of gravity offset are calculated using Fourier series expansion, representing the offset as the superposition of the fundamental frequency and its harmonics. The fundamental frequency is equal to the joint rotational speed. At the rated speed of 120 rpm (2 Hz), the amplitude of the fundamental frequency component of the center of gravity offset of joint 1 is... The value is 0.055 mm, the second harmonic component amplitude is 0.012 mm, and the third harmonic component amplitude is 0.0025 mm. The center of gravity offset is integrated over time using the trapezoidal integration method. The integration time is the work cycle time of 180 seconds, the integration step size is 0.01 seconds, and a total of 18,000 integration points are calculated. The instantaneous center of gravity offset is calculated at each integration point, and then the values ​​are accumulated to obtain the integral of the center of gravity offset. The integral of the center of gravity offset of joint 1 within the 180-second work cycle is 12.51 mm·s. The average integral of the center of gravity offset per unit time is calculated, which is equal to the total integral divided by the integration time. The average integral of the center of gravity offset of joint 1 is 12.51 / 180 = 0.0695 mm. The integral of the center of gravity offset is normalized by the maximum-minimum normalization method, with a normalization range of 0 to 1. The normalized integral of the center of gravity offset is 0.58. The data structure of the integral of the center of gravity offset is constructed, which includes information such as instantaneous center of gravity offset, radial and tangential components, time-varying characteristics, and integral.

[0077] Based on the calculated centrifugal vibration imbalance index and the obtained integral of center of gravity offset, the radial runout index during the operation of the robotic arm joint is derived. First, a radial runout physical model is established. Radial runout refers to the random displacement of the joint along the radial direction during operation, mainly caused by centrifugal vibration imbalance and center of gravity offset. The centrifugal vibration imbalance index of joint 1 is 0.139, and the normalized integral of center of gravity offset is 0.58. A radial runout index calculation model is established using a weighted summation method. The radial runout index is equal to the weighted sum of the centrifugal vibration imbalance index and the integral of center of gravity offset. The imbalance index weight is 0.4, and the weight of the centroid offset integral is 0.6. The radial runout index of joint 1 is calculated as 0.139×0.4+0.58×0.6=0.4036. A nonlinear mapping is performed on the radial runout index using an S-shaped function. The function parameters are set as follows: midpoint value 0.5, slope factor 8. The radial runout index is mapped to the range of 0 to 1, resulting in a mapped radial runout index of 0.42. The frequency characteristics of the radial runout are analyzed using power spectral density analysis to calculate the power spectrum of the radial runout. The main frequency components of the radial runout of joint 1 are... The power spectral densities for 2Hz (fundamental frequency), 4Hz (second harmonic), and 6Hz (third harmonic) are 0.025 mm² / Hz, 0.012 mm² / Hz, and 0.005 mm² / Hz, respectively. The statistical characteristics of radial runaway are calculated using probability density function analysis. The amplitude of the radial runaway follows a Rayleigh distribution with distribution parameters σ = 0.05 mm, a mean of 0.063 mm, a standard deviation of 0.033 mm, and a maximum value of 0.15 mm. Time series analysis of the radial runaway index is performed using an autoregressive moving average (ARMA) model. With parameters p=2 and q=1, the radial runaway index trend is predicted over the next 30 seconds. The prediction results show that the radial runaway index shows a slow upward trend, expected to reach 0.45 after 30 seconds. A radial runaway index data structure is constructed, including radial runaway index value, frequency characteristics, statistical characteristics, time series prediction results, etc. The data sampling rate is 100Hz, the time span is 180 seconds, and there are a total of 18,000 data points. Each data point includes a timestamp and a radial runaway index value. The radial runaway index data is output to provide basic data for subsequent force drift spatial interpolation processing.

[0078] Preferably, the spatial interpolation processing for the force drift of the robotic arm joints includes:

[0079] Time-varying characteristics of disordered triboacoustic wave intensity data are analyzed to output disordered time-varying intensity; energy distribution density skewness of disordered triboacoustic wave intensity data in different frequency bands is evaluated based on disordered time-varying intensity to obtain energy distribution density skewness data.

[0080] The radial traversal index is used to count the traversal frequency of the robotic arm joints, and then the difference in the amplitude of the radial offset is analyzed.

[0081] The acceleration change during the joint movement of the robotic arm is calculated based on the average difference of the displacement amplitude.

[0082] Based on the energy distribution density skewness data and the acceleration change, the spatial force offset of the robotic arm joint is quantified to obtain spatial force offset data.

[0083] Force drift spatial interpolation processing is performed on the spatial force offset data to obtain force drift interpolation data.

[0084] In this embodiment of the invention, time-varying characteristics of the obtained triboelectric acoustic wave disorder intensity data are analyzed. First, the disorder intensity data is rearranged according to the time series to form a time series dataset. The data sampling rate is 100Hz, the time span is 180 seconds, and there are a total of 18,000 data points. Each data point contains a timestamp and the corresponding disorder intensity value. The time series data is smoothed using a moving average filtering method with a window length of 5 points. The smoothness of the data is improved by 25% after filtering. Time-domain analysis is performed on the smoothed data to calculate the statistical characteristics of the disorder intensity, including mean 0.78, standard deviation 0.12, maximum value 0.92, minimum value 0.65, peak factor 3.25, skewness 0.45, and kurtosis 2.85. Time-frequency analysis is performed on the disorder intensity data using the Short Time Fourier Transform (STFT) method with a window length of 512 points, a window shift of 256 points, and a Hamming window as the window function. The time-frequency spectrum is calculated with a time resolution of 2.56 seconds and a frequency resolution of 0.195Hz. The main frequency components of the disordered intensity were extracted using a peak detection algorithm with a detection threshold set to 1.5 times the average amplitude. The detected main frequency components were 0.5Hz, 1.2Hz, 2.0Hz, and 3.5Hz, with corresponding amplitudes of 0.08, 0.06, 0.15, and 0.04, respectively. The periodicity characteristics of the disordered intensity were analyzed using autocorrelation analysis. The autocorrelation function was calculated with an autocorrelation length of 1000 points, corresponding to a time interval of 10 seconds. The autocorrelation function showed a significant peak at 0.5 seconds. The obvious peak value indicates that the disorder intensity has a short period of 0.5 seconds, and the secondary peak value at 2.0 seconds indicates that there is a long period of 2.0 seconds. Trend analysis of the disorder intensity data is carried out by using linear regression. The fitted trend line has a slope of 0.0005 / second, indicating that the disorder intensity generally shows a slight upward trend within 180 seconds, rising from 0.75 to 0.84. The disorder time-varying intensity data structure is constructed, which includes information such as time domain statistical characteristics, frequency domain characteristics, periodicity characteristics, and trend characteristics.

[0085] Based on the output disordered time-varying intensity data, the energy distribution density skewness of the disordered intensity data of friction sound waves in different frequency bands is evaluated. First, the frequency domain structure of the friction sound waves is divided into five frequency bands: low frequency (0-1kHz), mid-low frequency (1-3kHz), mid frequency (3-6kHz), mid-high frequency (6-10kHz), and high frequency (10-20kHz). Spectral data, including frequency points and corresponding amplitudes, is extracted for each band. The low frequency band contains 43 frequency points, the mid-low frequency band contains 86 frequency points, and the mid-high frequency band contains 86 frequency points. The frequency band contains 129 frequency points, the mid-high frequency band contains 172 frequency points, and the high frequency band contains 430 frequency points. The energy distribution of each frequency band is calculated using the power spectral density integral method to obtain the energy proportion of each band: low frequency band accounts for 5%, mid-low frequency band for 15%, mid frequency band for 45%, mid-high frequency band for 30%, and high frequency band for 5%. Statistical analysis is performed on the energy distribution of each frequency band, calculating the probability density function (PDF) of the energy distribution using the kernel density estimation method. A Gaussian kernel was selected as the function, and the bandwidth parameter was set to 0.1. PDF curves of energy distribution in each frequency band were obtained. The shape characteristics of the PDF curves were analyzed, and the skewness coefficients were calculated. Skewness is defined as the ratio of the third central moment to the cube of the standard deviation. The skewness was 0.25 for the low-frequency band, 0.42 for the mid-low-frequency band, 0.68 for the mid-frequency band, 0.55 for the mid-high-frequency band, and 0.32 for the high-frequency band. A correlation analysis was performed between the disordered time-varying intensity and the skewness of each frequency band. The Pearson correlation coefficient method was used to calculate the correlation between the disordered time-varying intensity and the skewness of each frequency band. The correlation coefficients with low-frequency skewness are 0.35, with mid-low-frequency skewness is 0.48, with mid-frequency skewness is 0.72, with mid-high-frequency skewness is 0.65, and with high-frequency skewness is 0.38. Based on the correlation analysis results, a weighted skewness model is constructed, with the weight of each frequency band proportional to the correlation coefficient. The weighted average skewness is calculated to be 0.56. An energy distribution density skewness data structure is constructed, which includes information such as the energy proportion of each frequency band, skewness value, correlation coefficient, and weighted average skewness.

[0086] Based on the derived radial runout index, the runout frequency of the robotic arm joints is counted. First, a runout frequency detection model is established, and the zero-crossing rate method is used to detect the runout frequency. The mean of the radial runout signal is subtracted to obtain a zero-mean signal. The number of times the signal crosses zero points per unit time is calculated. The runout frequency is equal to the number of zero crossings divided by 2. The sampling rate of the radial runout signal for joint 1 is 1000Hz, and the analysis window length is 1 second. A total of 180 windows are analyzed within a 180-second work cycle, with the runout frequency calculated once per window. In the first window, 36 zero crossings were detected, corresponding to a runout frequency of 18Hz. In the second window, 34 zero crossings were detected, corresponding to a runout frequency of... The frequency was calculated to be 17Hz, and so on. The average erratic frequency of 180 windows was 18.5Hz, with a standard deviation of 3.2Hz, a maximum value of 25Hz, and a minimum value of 12Hz. Spectral analysis of the erratic frequency was performed using Fast Fourier Transform (FFT). The FFT points were set to 1024, and the sampling rate was 1Hz, resulting in a frequency resolution of 0.001Hz. The main periods of change of the erratic frequency were 20 seconds, 60 seconds, and 120 seconds, with corresponding frequency components of 0.05Hz, 0.0167Hz, and 0.0083Hz. The radial displacement amplitude was analyzed using peak detection. The peak points of the surging signal were detected with a peak detection threshold set to the mean plus twice the standard deviation. A total of 3240 peak points were detected within 180 seconds, averaging 18 peak points per second, consistent with the surging frequency. The amplitude statistical characteristics of the peak points were calculated: the mean peak amplitude was 0.15 mm, the standard deviation was 0.06 mm, the maximum value was 0.28 mm, and the minimum value was 0.05 mm. The amplitude statistical characteristics of the trough values ​​were also calculated: the mean trough amplitude was -0.12 mm, the standard deviation was 0.05 mm, the maximum value was -0.04 mm, and the minimum value was -0.25 mm. The average displacement amplitude difference was calculated as the average of the absolute values ​​of the mean peak amplitude and the mean trough amplitude. The average displacement amplitude difference is calculated as (0.15 + 0.12) / 2 = 0.135 mm. A time series analysis of the average displacement amplitude difference is performed using a sliding window method with a window length of 10 seconds and a sliding step of 1 second. The average displacement amplitude difference within each window is calculated, and its trend is analyzed. The average displacement amplitude difference is relatively stable in the first 60 seconds, with an average value of 0.12 mm. It slightly increases between 60 and 120 seconds, with an average value of 0.135 mm. It increases significantly between 120 and 180 seconds, with an average value of 0.15 mm. A data structure for the turbulence frequency and average displacement amplitude difference is constructed, including information such as turbulence frequency, average displacement amplitude difference, and their time variation characteristics.

[0087] Based on the calculated average displacement amplitude difference, the acceleration change during the joint movement of the robotic arm is calculated. First, an acceleration calculation model is established. Based on Newton's second law, acceleration equals force divided by mass. The force during movement mainly comes from the impact force caused by the joint clearance. The impact force is proportional to the average displacement amplitude difference and the square of the movement frequency. The mass of joint 1 is 4.8 kg, the average displacement amplitude difference is 0.135 mm, and the movement frequency is 18.5 Hz. The average acceleration during movement is calculated using a simple harmonic motion model. The acceleration amplitude equals the displacement amplitude multiplied by the square of the angular frequency, and the angular frequency equals 2π multiplied by the movement frequency. The calculated acceleration amplitude is 1.83 m / s². Time-domain analysis of the acceleration is performed using a central difference method. Instantaneous acceleration was calculated using a differential step size of 0.001 seconds. The formula is: current acceleration equals the displacement at the next moment minus twice the current displacement plus the displacement at the previous moment, then divided by the square of the differential step size. A sequence of 180,000 instantaneous acceleration points over 180 seconds was calculated. The statistical characteristics of the acceleration were analyzed: mean acceleration was 0 m / s² (theoretically it should be 0, but the actual calculated result is 0.02 m / s², with the error mainly from numerical calculations); standard deviation was 0.75 m / s²; maximum value was 3.42 m / s²; minimum value was -3.25 m / s²; and peak factor was 4.56. Frequency domain analysis of the acceleration was performed using Fast Fourier Transform (FFT). For frequency components, the FFT points were set to 65536, the sampling rate was 1000Hz, and the frequency resolution was 0.015Hz. The main frequency components of acceleration were 18.5Hz (fundamental frequency), 37Hz (second harmonic), and 55.5Hz (third harmonic), with corresponding amplitudes of 1.5m / s², 0.6m / s², and 0.25m / s², respectively. The directional components of acceleration were calculated, decomposing them into x, y, and z components. Based on joint kinematics, the x-direction component accounted for 45%, the y-direction component for 40%, and the z-direction component for 15%. The mean x-direction acceleration was 0m / s², with a standard deviation of 0.34m / s². The mean y-direction acceleration was 0m / s². The mean acceleration in the z-direction is 0 m / s² with a standard deviation of 0.3 m / s². Time series analysis of acceleration variations was performed using a sliding window method with a window length of 1 second and a sliding step of 0.1 seconds. The root mean square (RMS) value of acceleration within each window was calculated, and its trend was analyzed. The RMS value was relatively stable in the first 60 seconds, with an average of 1.65 m / s². It slightly increased from 60 to 120 seconds, with an average of 1.85 m / s². It increased significantly from 120 to 180 seconds, with an average of 2.1 m / s². A data structure for acceleration variation was constructed, including information on acceleration time-domain characteristics, frequency-domain characteristics, directional components, and time-varying characteristics.

[0088] Based on the obtained energy distribution density skewness data and calculated acceleration changes, the spatial force offset of the robotic arm joint is quantified. First, the offset of the friction contact point is analyzed based on the energy distribution density skewness data. A regression analysis method is used to establish a mapping relationship between skewness and contact point offset. The mapping function uses a quadratic polynomial with coefficients [0.25, 0.35, 0.05]. The mid-frequency skewness is 0.68. The calculated friction contact point offset is 0.25 × 0.68² + 0.35 × 0.68 + 0.05 = 0.35 mm. The acceleration changes are then analyzed for directional components, extracting the acceleration components in the x, y, and z directions. The standard deviation of the acceleration in the x-direction is 0.34 m / s², and the standard deviation of the acceleration in the y-direction is... The acceleration is 0.3 m / s², the standard deviation of the acceleration in the z-direction is 0.11 m / s², and the calculated composite standard deviation of the acceleration is 0.47 m / s². The relationship between the acceleration direction component and the joint motion direction is analyzed. The main motion direction of joint 1 is rotation around the z-axis. The acceleration component in the xy-plane reflects the radial jaundice characteristics, and the acceleration component in the z-direction reflects the axial jaundice characteristics. The calculated radial jaundice acceleration is 0.45 m / s², and the axial jaundice acceleration is 0.11 m / s². Torsional torque offset is fitted based on the friction contact point offset and the jaundice acceleration direction component using the least squares method. The fitting model is a linear model. The torsional torque offset equals the friction contact point offset multiplied by the contact stiffness and then multiplied by the lever arm length. The contact stiffness is 2.5 × 10⁻⁶. 4 With a torque of N / m and a lever arm length of 0.035m, the calculated torsional moment offset is 0.35 × 10⁻⁶ N / m. -3 ×2.5×10 4 ×0.035=0.31N·m. Based on the torsional torque offset data and the axial acceleration direction component, a numerical simulation of the out-of-roundness of the minimum circumscribed circle of circular motion is performed. The least squares circle fitting algorithm is used, with 200 fitting points. The point coordinates are determined by the axial displacement and phase angle. The calculated minimum circumscribed circle radius is 0.25mm, the actual maximum radius of the trajectory is 0.28mm, the minimum radius is 0.22mm, and the out-of-roundness is 0.06mm. Based on the torsional torque offset data and the out-of-roundness simulation value, the spatial force offset of the robotic arm joint is quantified. The weighted average method is used, with the torsional torque offset weight being 0.6 and the out-of-roundness simulation value weight being 0.4. The calculated spatial force offset data has an average value of 0.31×0.6+0.06×0.4=0.21mm. A spatial force offset data structure is constructed, which includes information such as the friction contact point offset, the axial acceleration direction component, the torsional torque offset, the out-of-roundness simulation value, and the spatial force offset value.

[0089] Based on the obtained spatial force offset data, spatial interpolation processing of force drift is performed. First, the spatial force offset data is preprocessed, using median filtering to remove outliers. The filtering window length is 5 points, improving data smoothness by 20%. Spatial distribution analysis is then performed on the filtered data. The robotic arm's workspace is divided into a 10×10×10 grid, with each grid measuring 100mm×100mm×100mm. The force offset value at each grid point is analyzed. The average force offset in the central region of the workspace is 0.18mm, while the average force offset in the edge regions is 0.25mm. A spatial interpolation model is established, using radial basis function (RBF) interpolation. The interpolation method uses a quadratic function as the basis function, in the form f(r) = (r² + c²)^(1 / 2), where r is the spatial distance and c is the smoothing parameter, set to c = 0.5. 500 known points are selected as interpolation reference points, evenly distributed within the workspace. Each reference point contains spatial coordinates (x, y, z) and the corresponding force offset value. Interpolation calculations are performed on any point within the workspace, resulting in 5000 interpolation points, forming a 50×50×2 spatial grid with a grid spacing of 20 mm. The force offset value at each interpolation point is calculated using the formula: the force offset value at the interpolation point equals the force offset values ​​of all reference points multiplied by a factor of 1. The sum of corresponding weighting coefficients, determined by the radial basis function, is related to the distance between the interpolation point and the reference point; the closer the distance, the greater the weight. To evaluate the accuracy of the interpolation results, cross-validation is used. 100 points are randomly selected from 500 reference points as validation points, and the remaining 400 points are used for interpolation. The error between the interpolated value and the actual value is then calculated. The average absolute error is 0.015 mm, and the relative error is 7.1%. Spatial smoothing of the interpolated data is performed using a three-dimensional Gaussian filter with a kernel size of 3×3×3 and a standard deviation σ=1.0. The smoothed data is more continuous and reduces local fluctuations. A force drift interpolation data structure is then constructed. The data includes interpolation point coordinates, force offset interpolation values, and interpolation error estimates. It contains 5000 data points, each with spatial coordinates (x, y, z) and a corresponding force drift interpolation value. The output force drift interpolation data provides foundational data for subsequent nonlinear regression analysis. Furthermore, the interpolation data is visualized to generate a 3D force drift distribution cloud map. Colors from blue to red represent the increase in force drift value. The cloud map clearly shows the distribution pattern of force drift within the workspace, with smaller force drift in the central area and larger force drift in the edge areas, especially at the limit of the robotic arm's extension, where the force drift is greatest, reaching 0.32 mm.

[0090] Preferably, quantifying the spatial force offset of the robotic arm joints includes:

[0091] Analyze the offset of the friction contact point based on energy distribution density skewness data;

[0092] The directional component feature analysis of the acceleration change is performed to obtain the directional component of the surging acceleration;

[0093] Torsional torque offset is fitted based on the offset of friction contact point and the directional component of surging acceleration to obtain torsional torque offset data;

[0094] Based on the torsional torque offset data and the directional component of the surging acceleration, numerical simulation of the out-of-roundness of the minimum circumcircle of the circular motion is performed to obtain the out-of-roundness simulation value.

[0095] Based on the torsional torque offset data and out-of-roundness simulation values, the spatial force offset of the robotic arm joint is quantified to obtain spatial force offset data.

[0096] In this embodiment of the invention, based on the obtained energy distribution density skewness data, the offset of the friction contact point is analyzed. First, a physical model of skewness and contact point offset is established. The offset of the friction contact point leads to asymmetry in the energy distribution of the friction sound wave spectrum, manifested as an increase in spectral skewness. The skewness is 0.25 in the low-frequency band, 0.42 in the mid-low-frequency band, 0.68 in the mid-high-frequency band, 0.55 in the mid-high-frequency band, and 0.32 in the high-frequency band. A weighted average is calculated for the skewness of each frequency band, with the weight proportional to the energy proportion of each frequency band. The weight is 0.05 for the low-frequency band, 0.15 for the mid-low-frequency band, 0.45 for the mid-high-frequency band, 0.3 for the mid-high-frequency band, and 0.05 for the high-frequency band. The calculated weighted average skewness is 0.56. A mapping relationship between skewness and contact point offset is established. Using an experimental calibration method, the corresponding spectral skewness is measured under the condition of known contact point offset, obtaining multiple sets of skewness-offset data pairs. A skewness of 0.2 corresponds to... An offset of 0.1 mm corresponds to an offset of 0.22 mm for a skewness of 0.4, 0.38 mm for a skewness of 0.6, and 0.58 mm for a skewness of 0.8. Based on these data, a regression model is established using a quadratic polynomial fit, with fitting coefficients of [0.25, 0.35, 0.05] and a coefficient of determination R² of 0.95, indicating a good fit. Substituting the weighted average skewness of 0.56 into the regression model, the offset of the friction contact point is calculated to be 0.25 × 0.56² + 0.35 × 0.56 + 0.05 = 0.33 mm. Uncertainty analysis of the offset is performed using the Monte Carlo method, considering skewness measurement error and regression model error, and 10,000 simulations are conducted. The 95% confidence interval of the offset is obtained as [0.28 mm, 0.38 mm]. A data structure for the friction contact point offset is constructed, including the offset value, uncertainty estimate, regression model parameters, and other information.

[0097] To analyze the directional components of the calculated acceleration changes, a joint coordinate system is first established. The origin is located at the joint rotation center, the z-axis is along the joint rotation axis, the x-axis is along the joint radial direction, and the y-axis is determined according to the right-hand rule, forming a right-hand rectangular coordinate system. The acceleration vector is decomposed into the joint coordinate system, obtaining components in three directions. A coordinate transformation matrix is ​​used to transform the components. The transformation matrix is ​​determined by the current joint posture. When joint 1 is at 0°, the mean acceleration in the x-direction is 0 m / s², and the standard deviation is 0.34 m / s²; the mean acceleration in the y-direction is 0 m / s², and the standard deviation is 0.3 m / s²; and the mean acceleration in the z-direction is 0 m / s², and the standard deviation is 0.11 m / s². Spectral analysis is performed on the acceleration components in the three directions using Fast Fourier Transform (FFT). The FFT points are set to 4096 points, and the sampling rate is 1000 Hz, obtaining the frequency resolution. The acceleration frequency is 0.244 Hz, the dominant frequency of the acceleration in the x-direction is 18.5 Hz with an amplitude of 0.28 m / s², the dominant frequency of the acceleration in the y-direction is 18.5 Hz with an amplitude of 0.25 m / s², and the dominant frequency of the acceleration in the z-direction is 18.5 Hz with an amplitude of 0.09 m / s². The phase relationship of the acceleration components is analyzed, and the phase difference between the accelerations in the x and y directions is calculated. Using the cross-power spectral density method, the phase difference at the dominant frequency of 18.5 Hz is 90°, indicating that the acceleration in the xy plane exhibits circular motion characteristics. The radial composite acceleration is calculated, which is equal to the square root of the sum of the squares of the accelerations in the x and y directions. The mean of the radial acceleration is 0 m / s², the standard deviation is 0.45 m / s², and the maximum value is 1.85 m / s². A data structure for the directional components of the convective acceleration is constructed, containing information such as the statistical characteristics, spectral characteristics, phase relationship, and radial composite acceleration of the accelerations in the three directions. Based on the friction contact point offset obtained from step S7.1 and the directional component of the surging acceleration obtained from step S7.2, a torsional torque offset fitting is performed. First, a torsional torque offset model is established, which consists of two parts: friction torque offset and inertial torque offset. The friction torque offset is proportional to the friction contact point offset, with the proportionality coefficient being the contact stiffness multiplied by the lever arm length. The inertial torque offset is proportional to the surging acceleration, with the proportionality coefficient being the moment of inertia. The friction torque offset is calculated, with the friction contact point offset being 0.33 mm and the contact stiffness being 2.5 × 10⁻⁶ mm. 4 N / m, lever arm length is 0.035m, friction torque offset is 0.33×10 -3 ×2.5×10 4The inertial moment offset is calculated as follows: × 0.035 = 0.29 N·m. The standard deviation of radial acceleration is 0.45 m / s², and the moment of inertia of joint 1 is 0.08 kg·m². The standard deviation of the inertial moment offset is 0.08 × 0.45 = 0.036 N·m. The composite torsional moment offset is calculated using the mean square method. The standard deviation of the composite torsional moment offset is the square root of 0.29² + 0.036² = 0.292 N·m. Time series analysis of the torsional moment offset is performed using autoregression. A moving average ARMA model with order parameters p=2 and q=1 was used to analyze the time-varying characteristics of torsional torque offset. The torsional torque offset was relatively stable in the first 60 seconds with a mean of 0.28 N·m. It increased slightly from 60 to 120 seconds with a mean of 0.3 N·m, and increased significantly from 120 to 180 seconds with a mean of 0.35 N·m. A torsional torque offset data structure was constructed, which includes information such as friction torque offset, inertial torque offset, synthetic torsional torque offset, and their time-varying characteristics.

[0098] Based on the obtained torsional torque offset data and the analytically obtained axial acceleration direction component, a numerical simulation of the out-of-circle of the minimum circumcircle of the circular motion is performed. First, a circular motion model of the joint is established. Ideally, the joint rotation trajectory is a perfect circle. In reality, due to the influence of torsional torque offset and axial acceleration, the rotation trajectory deviates from the ideal circle, forming an irregular curve. Trajectory data for one rotation of the joint is collected at a sampling frequency of 1000Hz and a rotation speed of 60rpm, collecting 1000 points per rotation. Each point contains x and y coordinates. The x-coordinate is determined by the x-direction component of the axial displacement, and the y-coordinate... The displacement is determined by the y-direction component of the jaundice, which is obtained by quadratic integration of the jaundice acceleration using the trapezoidal rule with an integration step of 0.001 seconds. The mean displacement in the x-direction is 0 mm with a standard deviation of 0.12 mm, and the mean displacement in the y-direction is 0 mm with a standard deviation of 0.11 mm. The collected trajectory points are fitted with a minimum circumcircle using a least-squares circle fitting algorithm. The objective function is to minimize the sum of the squared distances from all points to the circle. The Levenberg-Marquardt algorithm is used for optimization, with 100 iterations and a convergence threshold of 10. -6The minimum circumcircle parameters were obtained by fitting, with the center coordinates at (0.02mm, 0.03mm) and a radius of 0.25mm. The distance from each trajectory point to the fitted circle was calculated, resulting in 1000 distance values. The mean distance was 0.05mm, the standard deviation was 0.03mm, the maximum distance was 0.12mm, and the minimum distance was 0.01mm. The out-of-roundness was calculated as the difference between the maximum and minimum distances, and was 0.11mm. Spectral analysis of the out-of-roundness was performed using a Fast Fourier Transform. Using FFT, the periodic characteristics of the out-of-roundness are analyzed. The FFT points are set to 1024, the sampling rate is 1000Hz, and the frequency resolution is 0.977Hz. The main frequency components of the out-of-roundness are 1Hz (fundamental frequency), 2Hz (second harmonic), and 3Hz (third harmonic), with corresponding amplitudes of 0.06mm, 0.03mm, and 0.01mm, respectively. A numerical data structure for the out-of-roundness simulation is constructed, including information such as the minimum circumcircle parameter, out-of-roundness, and spectral characteristics. The numerical values ​​of the out-of-roundness simulation are then output.

[0099] Based on the obtained torsional moment offset data and the simulated out-of-roundness values, the spatial force offset of the robotic arm joint is quantified. First, a spatial force offset model is established. The spatial force offset is jointly determined by the torsional moment offset and the out-of-roundness. The torsional moment offset causes uneven force distribution on the joint, while the out-of-roundness causes the joint motion trajectory to deviate from the ideal circle. Both factors jointly affect the spatial position accuracy of the joint. The angular deviation caused by the torsional moment offset is calculated; the angular deviation is equal to the torsional moment offset divided by the joint stiffness. The torsional moment offset is 0.292 N·m, and the joint stiffness is 2.8 × 10⁻⁶. 5 N·m / rad, angular deviation is 0.292 / 2.8×10 5 =1.04×10 -6 Given rad = 0.00006°, calculate the end effector position deviation caused by the angular deviation. The position deviation equals the angular deviation multiplied by the distance from the end effector to the joint. The distance from joint 1 to the end effector is 850 mm, and the position deviation is 850 × 1.04 × 10⁻⁶ mm. -6=0.000884mm, calculate the positional deviation caused by out-of-roundness. The positional deviation is equal to the out-of-roundness multiplied by the ratio of the distance from the end to the joint to the distance from the joint to the center of rotation. The out-of-roundness is 0.11mm, the distance from the joint to the center of rotation is 35mm, and the positional deviation is 0.11×850 / 35=2.67mm. Synthesize the spatial force offset using the weighted average method, with the torsional moment offset weighted at 0.3 and the out-of-roundness weighted at 0.7, and calculate the spatial force offset as 0.000884×0.3+2.67×0.7=1.87mm. Perform time series analysis on the spatial force offset using the sliding window method, with a window length of 1 second and a window sliding step of 0.1 seconds, and calculate each The mean spatial force offset within the window was analyzed, and its trend was examined. The spatial force offset was relatively stable in the first 60 seconds, with a mean of 1.8 mm. It increased slightly from 60 to 120 seconds, with a mean of 1.9 mm, and increased significantly from 120 to 180 seconds, with a mean of 2.1 mm. A spatial force offset data structure was constructed, including deviations caused by torsional moment offset, deviations caused by out-of-roundness, and information on the synthesized spatial force offset and its time-varying characteristics. The data sampling rate was 100 Hz, the time span was 180 seconds, and a total of 18,000 data points were collected. Each data point included a timestamp and the corresponding spatial force offset value. The spatial force offset data was output to provide basic data for subsequent spatial interpolation processing of force drift.

[0100] Preferably, step S3 includes:

[0101] The precision deviation quantification data is normalized to obtain the precision deviation normalized data.

[0102] Based on the accuracy deviation normalized data, the deviation increment is fitted in the time dimension, and the accuracy deviation increment fitted data is output.

[0103] Based on the accuracy deviation increment fitting data, adaptive training and correction of deviation are performed during the operation of the robotic arm to obtain adaptive correction data of deviation.

[0104] The random forest algorithm is used to design a deviation correction architecture for the adaptive deviation correction data, and the deviation correction architecture is sent to the terminal to execute the control of the robotic arm.

[0105] As an example of the present invention, reference is made to Figure 3 As shown, step S3 in this example includes:

[0106] S31: Normalize the precision deviation quantization data to obtain precision deviation normalized data;

[0107] In this embodiment of the invention, the obtained precision deviation quantification data is normalized. First, key features in the precision deviation quantification data are extracted, including joint angle deviation, end-effector position deviation, posture deviation, and comprehensive deviation index. The average angle deviation of joint 1 is 0.16°, the average angle deviation of joint 2 is 0.24°, the average angle deviation of joint 3 is 0.20°, the average x-axis component of end-effector position deviation is 0.85mm, the average y-axis component is 0.92mm, and the average z-axis component is 0.65mm, the average roll component of posture deviation is 0.28°, the average pitch component is 0.32°, and the average yaw component is 0.25°, and the comprehensive deviation index is 0.68. Maximum and minimum value normalization is performed on each feature. The normalization formula is: the normalized value equals the original value minus the minimum value, then divided by... The minimum and maximum values ​​are calculated by subtracting the maximum and minimum values. The minimum value for joint angle deviation is 0°, and the maximum value is set to 0.5°. The minimum value for end-effector position deviation is 0mm, and the maximum value is set to 2mm. The minimum value for attitude deviation is 0°, and the maximum value is set to 1°. The normalized joint angle deviation is calculated to be 0.32, joint angle deviation is 0.48, joint angle deviation is 0.40, end-effector position deviation x-axis component is 0.425, y-axis component is 0.46, z-axis component is 0.325, attitude deviation roll component is 0.28, pitch component is 0.32, yaw component is 0.25, and the normalized value of the comprehensive deviation index is 0.68. A precision deviation normalized data structure is constructed, which includes the normalized value and original value of each feature. The precision deviation normalized data is output to provide basic data for subsequent deviation increment fitting.

[0108] S32: Perform deviation increment fitting in the time dimension based on the precision deviation normalization data, and output the precision deviation increment fitting data;

[0109] In this embodiment of the invention, based on the obtained precision deviation normalized data, deviation increment fitting is performed in the time dimension. First, time series analysis is performed on the normalized data. The data sampling rate is 100Hz, the time span is 180 seconds, and a total of 18,000 data points are obtained. Each data point contains a timestamp and the corresponding precision deviation normalized value. The time series data is then downsampled with a sampling interval of 1 second, resulting in 180 data points after downsampling. Trend analysis is then performed on the downsampled data using the moving average method with a window length of 10 points. The moving average curve is calculated, and the trend of deviation over time is analyzed. The comprehensive deviation index rises from 0.58 to 0.65 in the first 60 seconds, from 0.65 to 0.72 in the 60 to 120 seconds, and from 0.7 in the 120 to 180 seconds. The deviation increment was increased from 2 to 0.82. A polynomial fitting method was used with a polynomial order of 3, a fitting time window of 180 seconds, and a time step of 1 second. The least squares method was used to determine the polynomial coefficients. The resulting incremental fitting function coefficients for the comprehensive deviation index were [0.0000012, 0.0000856, 0.0012, 0.58], with a determination coefficient R² of 0.96, indicating a good fitting effect. Similar fitting was performed on other features to obtain the incremental fitting function coefficients for joint angle deviation, end-effector position deviation, and posture deviation. A precision deviation incremental fitting data structure was constructed, containing the fitting function coefficients for each feature, fitting accuracy evaluation indicators, and other information. The precision deviation incremental fitting data was output, providing basic data for subsequent adaptive training and correction of deviation.

[0110] S33: Perform adaptive training and correction of deviations during the operation of the robotic arm based on the precision deviation increment fitting data to obtain adaptive correction data for deviations.

[0111] In this embodiment of the invention, based on the output precision deviation increment fitting data, adaptive training and correction of deviations during the robotic arm operation are performed. First, feature extraction is performed on the precision deviation increment fitting data using a sliding window method. The window length is 10 seconds, and the window sliding step is 2 seconds. Statistical features of the deviation within each window are calculated, including mean, variance, first-order difference, and second-order difference. The comprehensive deviation index has a mean of 0.59, a variance of 0.0004, a first-order difference mean of 0.0012, and a second-order difference mean of 0.0000025 within the first window (0-10 seconds). Based on the extracted features, a multi-dimensional feature vector is constructed with a dimension of 20, containing various deviation features and their statistics. Dimensionality reduction is performed on the feature vector using Principal Component Analysis (PCA), retaining principal components with a contribution rate of 95%. The dimension after dimensionality reduction is 8, resulting in a dimensionality-reduced feature representation. A gated recurrent unit for deviation correction is then performed based on this dimensionality-reduced feature representation. The time-series iterative training of the GRU network was performed. The network structure included 8 input layer nodes, 32 hidden layer nodes, and 6 output layer nodes. The activation function was tanh, the initial bias of the forget gate was set to 1.0, the learning rate was 0.001, the batch size was 16, the training epochs were 200, the training dataset size was 150 samples, and the validation set size was 30 samples. An early stopping strategy was adopted during training; training was stopped if the loss did not decrease after 10 consecutive validation epochs to obtain corrected and optimized training data. The root mean square error of the model on the validation set was 0.025. Based on the trained GRU model, adaptive training and correction of deviations during the operation of the robotic arm were performed. The current state features were input in real time, and the model output position corrections of x-axis -0.42mm, y-axis -0.46mm, z-axis -0.32mm, and angle corrections of roll -0.14°, pitch -0.16°, and yaw -0.12° to obtain the adaptively corrected deviation data.

[0112] S34: Use the random forest algorithm to design a deviation correction architecture for the adaptive deviation correction data, and send the deviation correction architecture to the terminal to execute the control of the robotic arm.

[0113] In this embodiment of the invention, based on the bias adaptive correction data obtained in step S3.3, a bias correction architecture is designed using the random forest algorithm. First, the correction data is preprocessed: feature values ​​are standardized, outliers are removed, and the dataset is divided into training and test sets, with the training set accounting for 80% and the test set for 20%. A random forest model is constructed, with the following parameters set: 100 trees, a maximum depth of 15, a minimum number of leaf node samples of 5, and the feature splitting criterion being the root mean square error. The number of randomly selected features is the square root of the total number of features. The random forest model is trained using the training set, and 5-fold cross-validation is used to evaluate the model's performance during training. The model's root mean square error on the test set is 0.018, and the coefficient of determination (R²) is 0.97, indicating high prediction accuracy and effective feature analysis. The importance of positional deviation features is 0.45, the importance of angle deviation features is 0.35, and the importance of time features is 0.2. A deviation correction architecture is constructed based on a random forest model, which includes a feature extraction module, a feature dimensionality reduction module, a prediction module, and an execution module. The correction architecture is compiled into a binary file with a size of 2.5MB and sent to the terminal controller via TCP / IP protocol at a transmission rate of 10Mbps. The terminal controller model is TC-500, with a processor clock frequency of 1.2GHz and 2GB of memory. After receiving the correction architecture, the terminal loads it into memory, starts the execution thread, and the sampling frequency is 100Hz. The correction command latency is less than 5ms, which executes precise control of the robotic arm, thereby improving positional accuracy and trajectory smoothness.

[0114] Preferably, the adaptive training correction of deviations during the robotic arm operation includes:

[0115] Feature extraction is performed on the precision deviation increment fitting data, and the statistical characteristics of the deviation, including mean and variance, are calculated using the sliding window method to obtain the deviation feature data.

[0116] Based on the deviation feature data, a multidimensional feature vector is constructed, and the feature vector is then subjected to dimensionality reduction processing to obtain the dimensionality-reduced feature representation;

[0117] Time-series iterative training of a gated recurrent unit network based on the dimensionality-reduced feature representation is used to correct biases and obtain corrected and optimized training data.

[0118] Based on the corrected and optimized training data, adaptive training correction of deviations during the operation of the robotic arm is performed to obtain adaptive correction data for deviations.

[0119] In this embodiment of the invention, feature extraction is performed on the output precision deviation increment fitting data. A sliding window method is used to calculate the statistical characteristics of the deviation. The window length is set to 10 seconds, and the window sliding step size is set to 2 seconds, forming 86 windows in 180 seconds of data. Statistical characteristics are calculated for the data within each window, including mean, variance, maximum value, minimum value, peak factor, skewness, kurtosis, first-order difference mean, and second-order difference mean. Taking the first window (0-10 seconds) as an example, the comprehensive deviation index has a mean of 0.59, a variance of 0.0004, a maximum value of 0.61, a minimum value of 0.58, a peak factor of 3.2, a skewness of 0.35, a kurtosis of 2.8, a first-order difference mean of 0.0012, and a second-order difference mean of 0.0000025. The same statistical characteristics are calculated for joint angle deviation, end-effector position deviation, and attitude deviation. The mean of joint 1 angle deviation within the first window is 0.28, and the variance is... The variance is 0.0002. The mean and variance of the angle deviation of joint 2 are 0.42, and the mean and variance of the angle deviation of joint 3 are 0.35, and the variance of the angle deviation of joint 3 are 0.00025. The mean and variance of the x-axis component of the end-effector position deviation are 0.38, and the variance of the x-axis component are 0.0005, the mean and variance of the y-axis component are 0.41, and the variance of the y-axis component are 0.0006, the mean and variance of the z-axis component are 0.29, and the variance of the z-axis component are 0.0004. The mean and variance of the roll component of the attitude deviation are 0.25, and the variance of the roll component are 0.0003, the mean and variance of the pitch component are 0.28, and the variance of the pitch component are 0.00035. The mean and variance of the yaw component are 0.22, and the variance of the pitch component are 0.00025. The statistical features of all windows are combined to form deviation feature data. The data dimension is 86×45, that is, 86 windows, each window has 45 features. A deviation feature data structure is constructed, which includes the statistical feature values ​​of each window. The deviation feature data is output, providing basic data for the subsequent construction of multi-dimensional feature vectors.

[0120] Based on the obtained deviation feature data, a multidimensional feature vector is constructed. First, the feature data is standardized using the Z-score standardization method, subtracting the mean of each feature and dividing by its standard deviation, resulting in a mean of 0 and a standard deviation of 1 for all features. The standardized feature values ​​range approximately from -3 to 3. Correlation analysis is then performed on the standardized features, calculating the Pearson correlation coefficient between features. Highly correlated features with an absolute correlation coefficient greater than 0.9 are removed, retaining more representative features. The number of features is reduced from 45 to 32. A multidimensional feature vector with dimension 32 is then constructed, containing statistical features of comprehensive deviation indices, joint angle deviations, end-effector position deviations, and posture deviations. Dimensionality reduction is then performed on the multidimensional feature vector using principal component analysis (PCA). First, the feature covariance matrix is ​​calculated. The original feature vectors are transformed using a 32×32 matrix. The eigenvalues ​​and eigenvectors of the covariance matrix are calculated, with the eigenvalues ​​sorted from largest to smallest as follows: 5.8, 4.2, 3.1, 2.5, 1.9, 1.4, 1.1, 0.8, 0.6, 0.5, etc. The contribution rate and cumulative contribution rate of each principal component are calculated. The cumulative contribution rate of the first 8 principal components is 95.2%, exceeding the preset 95% threshold. Therefore, the first 8 principal components are retained. A PCA transformation matrix of size 32×8 is constructed, and a linear transformation is performed on the original feature vectors to obtain the dimensionality-reduced feature representation. The feature dimension after dimensionality reduction is 8, reducing the data volume by 75%. A data structure for the dimensionality-reduced feature representation is constructed, containing the dimensionality-reduced eigenvalues, the PCA transformation matrix, the eigenvalue contribution rate, etc. The dimensionality-reduced feature representation is output, providing basic data for subsequent training of the gated recurrent unit network.

[0121] Based on the obtained dimensionality-reduced feature representations, a gated recurrent unit (GRU) network for bias correction is subjected to time-series iterative training. First, a training dataset is constructed. The input consists of the feature representations from five consecutive time steps, and the output is the bias correction for the next time step, including position and angle corrections, for a total of six output values. The dataset size is 80 samples, divided into a training set of 64 samples and a validation set of 16 samples. The GRU network structure is designed with 8 nodes in the input layer (corresponding to the dimensionality of the reduced features), a single-layer GRU structure with 32 nodes in the hidden layer, and 6 nodes in the output layer (corresponding to the six correction values). The total number of network parameters is 4518. Training parameters are set, the loss function is mean squared error (MSE), the optimization algorithm is Adam, the initial learning rate is 0.001, and a learning rate decay strategy is adopted, decreasing to 0.8 times the original value every 50 epochs. The batch size is 16, and the maximum number of training epochs is set to 500. During network training, after the first round of training, the training set loss was 0.185 and the validation set loss was 0.192. After the 50th round of training, the training set loss decreased to 0.068 and the validation set loss decreased to 0.075. After the 100th round of training, the training set loss decreased to 0.042 and the validation set loss decreased to 0.048. After the 150th round of training, the training set loss decreased to 0.031 and the validation set loss decreased to 0.035. After the 200th round of training, the training set loss decreased to 0.025 and the validation set loss decreased to 0.028. After that, the validation set loss no longer decreased significantly, triggering the early stopping strategy, stopping training, and saving the trained model parameters. The root mean square error of the model on the validation set was 0.167 and the mean absolute error was 0.132. A corrected and optimized training data structure was constructed, including model parameters, training history, performance evaluation metrics, and other information. The corrected and optimized training data was output to provide basic data for subsequent bias adaptive training correction.

[0122] Based on the obtained corrected and optimized training data, adaptive training correction of deviations during the robotic arm operation is performed. First, the trained GRU model parameters are loaded to construct a real-time prediction system. The system sampling frequency is 100Hz, collecting robotic arm state data every 10ms, including joint angles, end effector position, and posture. Feature extraction is performed on the collected state data using the same sliding window method as in the training phase, with a window length of 10 seconds and a sliding step of 0.1 seconds. Statistical features are calculated, and the extracted features are standardized using Z-score standardization with the mean and standard deviation saved during the training phase. Dimensionality reduction is then performed on the standardized features using the PCA transformation matrix saved during the training phase, reducing the 32-dimensional features to 8-dimensional features. The dimensionality-reduced features are input into the GRU model, and the model outputs position and angle corrections. The position correction is shown on the x-axis. The model output corrections are -0.42mm (y-axis), -0.46mm (z-axis), and -0.32mm (roll, pitch, yaw). Angle corrections of -0.14°, -0.16°, and -0.12° are applied. An exponentially weighted moving average method with a smoothing coefficient α=0.3 is used to reduce abrupt changes in the corrections. The smoothed corrections are then applied to the robotic arm control system to correct the robotic arm's motion trajectory, achieving adaptive compensation for accuracy deviations. Accuracy data before and after correction are recorded. Before correction, the average end-effector position error was 0.85mm, which decreased to 0.25mm after correction. The average attitude error was 0.28°, which decreased to 0.08° after correction. An adaptive deviation correction data structure is constructed, containing information such as correction amounts and correction effect evaluation indicators. This outputs adaptive deviation correction data, providing foundational data for subsequent deviation correction architecture design.

[0123] The present invention also provides a control system for a robotic arm, for executing the control method for the robotic arm described above, the control system comprising:

[0124] The friction acoustic wave recognition module is used to deploy miniature acoustic wave sensors at the joints of a robotic arm and monitor the acoustic waves during the operation of the robotic arm joints to obtain the joint operation monitoring acoustic waves; the joint operation monitoring acoustic waves are marked with friction acoustic wave recognition to output the frequency domain structure of friction acoustic waves;

[0125] The precision deviation quantization module is used to derive the radial displacement index of the robotic arm joint during operation based on the frequency domain structure simulation of the friction acoustic wave, and then perform spatial interpolation processing of the force drift of the robotic arm joint to output force drift regression data; based on the force drift regression data, the running precision deviation of the joint is quantified to obtain precision deviation quantization data.

[0126] The correction architecture design module is used to perform adaptive training and correction of deviations during the operation of the robotic arm based on the precision deviation quantification data, obtain deviation adaptive correction data, and then design the deviation correction architecture to build the deviation correction architecture. The deviation correction architecture is then sent to the terminal to execute the control of the robotic arm.

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

Claims

1. A control method for a robotic arm, characterized in that, Includes the following steps: Step S1: Deploy miniature acoustic sensors at the joints of the robotic arm and monitor the acoustic waves during the operation of the robotic arm joints to obtain the joint operation monitoring acoustic waves; Friction acoustic wave identification and labeling are performed on the joint operation monitoring acoustic waves to output the frequency domain structure of friction acoustic waves; Step S2: Based on the frequency domain structure simulation of the friction sound wave, derive the radial displacement index of the robotic arm joint during operation, and then perform spatial interpolation processing of the force drift of the robotic arm joint to output force drift regression data; quantify the running accuracy deviation of the joint based on the force drift regression data to obtain accuracy deviation quantification data. Step S3: Based on the precision deviation quantification data, perform deviation adaptive training and correction during the operation of the robotic arm to obtain deviation adaptive correction data. Then, design the deviation correction architecture to build the deviation correction architecture and send the deviation correction architecture to the terminal to execute the control of the robotic arm. Step S2 includes: The frequency domain structure of friction acoustic waves is quantized to generate disordered intensity data of friction acoustic waves. The radial runout index during the operation of the robotic arm joint is derived by simulating the random intensity data of frictional acoustic waves. Based on the disordered intensity data of frictional acoustic waves and the radial runout index, the force drift space interpolation of the robotic arm joint is processed to obtain the force drift interpolation data. Nonlinear regression analysis is performed on the force drift interpolation data to output force drift regression data; Based on the force drift regression data and radial runout index, the running accuracy deviation of the joint is quantified to obtain the accuracy deviation quantification data. Disorder intensity quantization includes: The narrowband frequency sideband amplitude ratio of the triboacoustic wave frequency domain structure is calculated to obtain the frequency narrowband sideband amplitude ratio. Based on the frequency narrowband sideband amplitude ratio, frequency amplitude abrupt fractal features of the frequency domain structure of the friction acoustic wave are identified, and amplitude abrupt fractal features are obtained. The peak-valley distribution density interval is calculated for the fractal characteristics of amplitude abrupt change, and then the height ratio variance of the peaks is calculated. Based on the high proportional variance, the disordered entropy value of the amplitude abrupt fractal feature is differentiated by the disordered entropy value to obtain the frequency disordered entropy value. The disorder intensity is quantized based on the frequency disorder entropy value and the frequency narrowband sideband amplitude ratio to generate triboacoustic wave disorder intensity data. The radial runout index during the operation of the robotic arm joints is derived as follows: Obtain the structural design, mass, and stiffness of the robotic arm joints; obtain the joint clearance at the time of manufacture based on the structural design of the robotic arm joints; Based on the disordered intensity data of frictional acoustic waves, the mass and stiffness of the joint are simulated by resonance intensity coupling to obtain the joint resonance coupling intensity. Based on the aforementioned friction acoustic wave disorder intensity data and joint resonance coupling intensity, a random process simulation of joint gap reciprocating expansion at the factory is performed, and joint gap expansion data is output. The eccentric vibration imbalance index of joint motion is determined based on joint space enlargement data and joint resonance coupling strength. The center of gravity offset integral is obtained by integrating the center of gravity offset of the joint movement based on the centrifugal vibration imbalance index; the radial movement index of the robotic arm joint during operation is derived based on the centrifugal vibration imbalance index and the center of gravity offset integral. Spatial interpolation processing for force drift of robotic arm joints includes: Time-varying characteristics of disordered triboacoustic wave intensity data are analyzed to output disordered time-varying intensity; energy distribution density skewness of disordered triboacoustic wave intensity data in different frequency bands is evaluated based on disordered time-varying intensity to obtain energy distribution density skewness data. The radial traversal index is used to count the traversal frequency of the robotic arm joints, and then the difference in the amplitude of the radial offset is analyzed. The acceleration change during the joint movement of the robotic arm is calculated based on the average difference of the displacement amplitude. Based on the energy distribution density skewness data and the acceleration change, the spatial force offset of the robotic arm joint is quantified to obtain spatial force offset data. Force drift spatial interpolation processing is performed on the spatial force offset data to obtain force drift interpolation data; The spatial force offset quantization of the robotic arm joints includes: Analyze the offset of the friction contact point based on energy distribution density skewness data; The directional component feature analysis of the acceleration change is performed to obtain the directional component of the surging acceleration; Torsional torque offset is fitted based on the offset of friction contact point and the directional component of surging acceleration to obtain torsional torque offset data; Based on the torsional torque offset data and the directional component of the surging acceleration, numerical simulation of the out-of-roundness of the minimum circumcircle of the circular motion is performed to obtain the out-of-roundness simulation value. Based on the torsional torque offset data and out-of-roundness simulation values, the spatial force offset of the robotic arm joint is quantified to obtain spatial force offset data.

2. The control method for the robotic arm according to claim 1, characterized in that, Step S1 includes the following steps: Miniature acoustic sensors are deployed at the joints of the robotic arm to monitor the acoustic waves during the joint operation, thus obtaining the acoustic waves for joint operation monitoring. The acoustic waves used for joint movement monitoring are processed into frames to generate monitoring frame acoustic waves. Frictional sound waves are obtained by identifying and marking the monitored frame-by-frame sound waves. The friction sound wave is frequency-domain converted to output the frequency domain structure of the friction sound wave.

3. The control method for the robotic arm according to claim 1, characterized in that, Step S3 includes: The precision deviation quantification data is normalized to obtain the precision deviation normalized data. Based on the accuracy deviation normalized data, the deviation increment is fitted in the time dimension, and the accuracy deviation increment fitted data is output. Based on the accuracy deviation increment fitting data, adaptive training and correction of deviation are performed during the operation of the robotic arm to obtain adaptive correction data of deviation. The random forest algorithm is used to design a deviation correction architecture for the adaptive deviation correction data, and the deviation correction architecture is sent to the terminal to execute the control of the robotic arm.

4. The control method for the robotic arm according to claim 3, characterized in that, Adaptive training and correction of deviations during robotic arm operations include: Feature extraction is performed on the precision deviation increment fitting data, and the statistical characteristics of the deviation, including mean and variance, are calculated using the sliding window method to obtain the deviation feature data. Based on the deviation feature data, a multidimensional feature vector is constructed, and the feature vector is then subjected to dimensionality reduction processing to obtain the dimensionality-reduced feature representation; Time-series iterative training of a gated recurrent unit network based on the dimensionality-reduced feature representation is used to correct biases and obtain corrected and optimized training data. Based on the corrected and optimized training data, adaptive training correction of deviations during the operation of the robotic arm is performed to obtain adaptive correction data for deviations.

5. A control system for a robotic arm, characterized in that, For performing the control method of the robotic arm as described in claim 1, the control system of the robotic arm includes: The friction acoustic wave recognition module is used to deploy miniature acoustic wave sensors at the joints of a robotic arm and monitor the acoustic waves during the operation of the robotic arm joints to obtain the joint operation monitoring acoustic waves; the joint operation monitoring acoustic waves are marked with friction acoustic wave recognition to output the frequency domain structure of friction acoustic waves; The precision deviation quantization module is used to derive the radial displacement index of the robotic arm joint during operation based on the frequency domain structure simulation of the friction acoustic wave, and then perform spatial interpolation processing of the force drift of the robotic arm joint to output force drift regression data; based on the force drift regression data, the running precision deviation of the joint is quantified to obtain precision deviation quantization data. The correction architecture design module is used to perform adaptive training and correction of deviations during the operation of the robotic arm based on the precision deviation quantification data, obtain deviation adaptive correction data, and then design the deviation correction architecture to build the deviation correction architecture. The deviation correction architecture is then sent to the terminal to execute the control of the robotic arm.

Images