The robot teleoperation system is based on two IMUs and an operating handle. The operating handle provides the position and orientation of the human hand, while the two IMUs, worn on the upper arm and forearm of the human body, are used to detect the motion state of the human arm. The data from the IMUs and operating handle are transmitted via Bluetooth to the STM32 for data acquisition and processing. Subsequently, the data are sent through a Wi-Fi module to the rk3588 processor for robot inverse kinematics and arm angle calculations, generating the trajectory of a seven-degree-of-freedom (7-DOF) robot. The system composition and data transmission process are illustrated in Figure 1.

This paper adopts incremental position mapping and absolute attitude mapping. The kinematic model of the robot and the human body is established before determining the human–robot motion mapping. The 7-DOF robot structure is shown in Figure 2 and the robot DH parameters are provided in Table 1.

The robot kinematics model can be obtained by: T r θ 1 , θ 2 , θ 3 , θ 4 , θ 5 , θ 6 , θ 7 7 0 = r 11 r 12 r 13 p x r 21 r 22 r 23 p y r 31 0 r 32 0 r 33 p z 0 1 = R r 7 0 P 7 0 0 1 . (1)

Here, P 7 0 = p r , x , p r , y , p r , z is the robot end-effector position.

To solve the problem of multiple solutions in the inverse solution of a redundant DOF manipulator, a unique solution can be obtained by introducing the global configuration parameter GCk. According to Faria et al. (2018), the spatial orientation of robotic arms is directly affected by GCk of the shoulder, elbow, and wrist.

To address the global and local self-motion manifolds, two supplementary parameters are incorporated into the redundant robot inverse kinematics calculation, namely, global configuration (GC) and arm angle ψ . The global configuration is used to specify the branch of the inverse kinematics solutions for the global configuration manifold. The arm angle ψ indicates the elbow position in the redundancy circle, as shown in Figure 3. The global configuration parameter G C k is divided into three variables that represent the sign of the joint angle coordinates for the shoulder joint ( G C 2 ), the elbow joint ( G C 4 ), and the wrist joint ( G C 6 ). G C k is given as follows: G C k = 1 , if θ k ≥ 0 − 1 , if θ k < 0 . , ∀ k ∈ 2 , 4 , 6 . (2)

In this paper, when the operator employs the right arm for teleoperation, the values are set as GC2 = 1, GC4 = 1, and GC6 = −1; conversely, when the left arm is utilized for teleoperation, the values are assigned as GC2 = −1, GC4 = −1, and GC6 = 1. Subsequently, a unique inverse solution can be derived based on the pose T s , t 1 , the arm angle ψ , and the global configuration parameter GCk. Since θ_4 is not affected by the arm angle ψ , other angles can be obtained using the following formulas. θ 1 = atan ⁡ 2 G C 2 a s 22 ⁡ sin ⁡ ψ + b s 22 ⁡ cos ⁡ ψ + c s 22 , G C 2 a s 12 ⁡ sin ⁡ ψ + b s 12 ⁡ cos ⁡ ψ + c s 12 , θ 2 = G C 2 ⁡ arccos ⁡ a s 32 ⁡ sin ⁡ ψ + b s 32 ⁡ cos ⁡ ψ + c s 32 , θ 3 = atan ⁡ 2 G C 2 − a s 33 ⁡ sin ⁡ ψ − b s 33 ⁡ cos ⁡ ψ − c s 33 , G C 2 − a s 31 ⁡ sin ⁡ ψ − b s 31 ⁡ cos ⁡ ψ − c s 31 , θ 5 = atan ⁡ 2 G C 6 a w 23 ⁡ sin ⁡ ψ + b w 23 ⁡ cos ⁡ ψ + c w 23 , G C 6 a w 13 ⁡ sin ⁡ ψ + b w 13 ⁡ cos ⁡ ψ + c w 13 , θ 6 = G C 6 ⁡ arccos ⁡ a w 33 ⁡ sin ⁡ ψ + b w 33 ⁡ cos ⁡ ψ + c w 33 , θ 7 = atan ⁡ 2 G C 6 a w 32 ⁡ sin ⁡ ψ + b w 32 ⁡ cos ⁡ ψ + c w 32 , G C 6 − a w 31 ⁡ sin ⁡ ψ − b w 31 ⁡ cos ⁡ ψ − c w 31 . (3)

To construct a virtual human arm model, the coordinate system on the human arm is established as shown in Figure 1. Coordinate systems {1} and {2} are on the human shoulder joint, while coordinate system {3} is on the human elbow joint, and coordinate system {4} is established on the wrist. The DH parameterof virtual huma arm is shown in Table 2.

Here, a h , 2 is the length of the upper arm and a h , 3 is the length of the lower arm. According to the homogeneous transformation formula, the elbow position can be derived as follows: P h 3 0 = a h , 2 c h , 2 c h , 1 a h , 2 c h , 2 s h , 1 a h , 2 s h , 2 . (4)

The wrist position is P h 4 0 = c h , 1 a h , 3 c h , 23 + a h , 2 c h , 2 s h , 1 a h , 3 c h , 23 + a h , 2 c h , 2 a h , 3 s h , 23 + a h , 2 s h , 2 . (5)

This paper employs a hybrid mapping approach, where position is mapped incrementally while orientation uses absolute mapping. The robot model and virtual human arm model are utilized to map the posture of both the human arm and wrist.

The positions of human wrist and robot can be obtained according to Equations 1, 5. The positions of the human wrist and robot end-effector are mapped incrementally. Assuming that the initial calibration position is P m , t 0 = P h , t 0 and P s , t 0 = P r , t 0 7 0 , the positions in relation to the master and slave ends can be expressed as follows: P h , t 1 = P h , t 0 + ∆ P h , t 1 , P r , t 1 = P r , t 0 + ∆ P r , t 1 , (6) ∆ P r , t 1 = K ∆ P h , t 1 .

Here, K is the mapping scale. The attitude of the human wrist and the robot end-effector are mapped absolutely: ∆ R r , t 1 = ∆ R h , t 1 . (7)

In order to ensure that the robot arm can adjust the attitude of the intermediate joint according to the requirements of the operator to achieve obstacle avoidance, the intermediate joint adopts absolute attitude mapping: T h , t 1 = P h , t 1 , ∆ R h , t 1 , T r , t 1 = K ∆ P r , t 1 , ∆ R r , t 1 . (8)

According to Equations 6–8, the relationship mapping can be established. For ease of calculation, the arm angle ψ in this paper is defined as the angle between the arm and reference planes, where the arm plane is considered parallel to the axis X 1 and X 2 of the IMU on the upper and lower arms. To figure out the arm plane, the Euler angles of the two IMUs are converted into direction vectors. By defining the IMU Euler angles as α 1 , β 1 , γ 1 and α 2 , β 2 , γ 2 , the normal vector S of the arm plane can be obtained by the cross product of X 1 and X 2 . S = X 1 × X 2 = c γ 1 * c β 1 s γ 1 * c β 1 s β 1 × c γ 2 * c β 2 s γ 2 * c β 2 s β 2 . (9)

Defining the reference plane as being perpendicular to the horizontal plane, the arm angle ψ can be obtained by dotting the normal vector S with the horizontal normal vector H . ψ = a r c o s S · H . (10)

The inverse kinematics can be obtained by substituting Equation 10 into Equation 2.

When there is magnetic field interference in the operating environment, the magnetometer cannot be used to calibrate the IMU online, and the IMU cumulative error e 1 t , e 2 t will appear. In this case, the normal vector of the arm plane is S ′ e 1 t , e 2 t , and the arm angle error is ψ ′ t = a r ⁡ cos S ′ t · H . (11)

To observe the drift problem, an IMU with magnetometer malfunction is bound with a high-precision IMU with magnetometer. A random arm motion is conducted to observe the drift of the low-cost IMU, while the high-precision IMU is regarded as a reference standard. The Figure 4 shows the IMU yaw angle without (blue line) and with the magnetometer after 20 min of various random motions. The final cumulative error is approximately 70°. It is worth noting that to perform online calibration of the IMU on the arm, accurate data are required as a reference point. Since the handle has its own vision-based calibration function, this paper assumes that the position data of the handle are accurate.

Considering the computation and time costs for online training, random forest-based feature selection is conducted to exclude unimportant feature data. Random forest-based feature selection combines the powerful representation capabilities of neural networks with the stability and nonlinear processing capabilities of random forest models. During the generation of each decision tree, the data subset and feature subset are randomly selected to enhance the robustness of the model and reduce overfitting. The input data for feature selection include the position and attitude of the operating handle, the attitude of IMUs, and the angular velocity and acceleration of IMUs. The feature selection results indicate that the position of the operating handle and the pitch and roll angle of IMUs have a significant impact on the IMUs’ yaw angle estimation.

In this paper, Kalman filtering is used to reduce IMU data bias and noise. To achieve this, the angular velocity measured by the gyroscope is used to predict the attitude and simultaneously model the slow change of zero offset: x k = A x k − 1 + B u k + w k . (12)

Here, A is the state transition matrix, u k is the angular velocity measured by the gyroscope, and w k is data noise. The process is started with the initiating states x 0 and P 0 . The prediction state x ^ k − and error covariance P k − are calculated as follows: x ^ k − = A x ^ k − 1 + B u k , P k − = A P k − 1 A ∼ T + Q . (13)

The Kalman gain K k and state are updated according to K k = P k − H T H P k − H T + R − 1 , x ^ k = x ^ k − + K k Z k − H x ^ k − , P k = P k − − K k H P k − . (14)

Here, Q is the process noise covariance matrix, and R is the measurement noise covariance matrix.

The ML-LSTM neural network introduces the modulation gate into traditional LSTM to evaluate the importance of historical information at different times. It can improve the traditional memory mode and take the memory summation average value of each data segment as the standard memory, thus improving the model’s ability to recognize key information and the model’s ability to interpret. The structure of the ML-LSTM neural network is shown in Figure 6.

The forget gate of ML-LSTM is used to determine what information can pass through the memory unit and generate an f k according to the output value h k − 1 at the last moment and the current input value x k . f k = σ W f · h k − 1 , x k + b f . (15)

Here, W f is the weight, b f is the offset, and σ is the activation function sigmoid. To generate updated information, the values of i k and the new candidate C ∼ k will be calculated and added to the memory unit as a candidate value generated by the current layer. i k = σ W i · h k − 1 , x k + b i , C ∼ k = tanh ⁡ W C · h k − 1 , x k + b C . (16)

The memory cells are updated by C k ′ = f k · C k − 1 + i k · C ∼ k . (17)

The modulation gate sums the memory information and calculates the mean E. C k = E = ∑ 1 m   C k , m ′ m . (18) Here, m is the number of data segments. The model output h k can be obtained as follows: h k = σ W o · h k − 1 , x k + b o · ⁡ tanh C k . (19)

Then, the PSO algorithm is run offline using previous data to optimize the neuron number n n n and learning rate l r of ML-LSTM. The initial position and speed of the particles are first randomly generated. The fitness evaluation of each particle solution is conducted. The historical best position ( p i d t ) and the global best position ( p g d t ) for each particle are then recorded. The speed and position of each particle according to p i d t and p g d t are updated to calculate the neuron number and learning rate according to the following equation. v i d t + 1 = ω v i d t 1 + c 1 r 1 p i d t − x i d t + c 2 r 2 p g d t − x g d t , x i d t + 1 = x i d t + v i d t + 1 . (20)

Here, v i d t + 1 and x i d t + 1 can be represented by v l r and v n n and l r and n n n . The above process will be repeated until the number of iterations is reached or an optimal solution is found. The PSO flow chart is shown in Figure 7.

To evaluate the performance of the IMU state estimation model, three evaluation indexes are used, including the root mean square error (RMSE), mean absolute error (MAE), mean bias error (MBE), and R-squared R 2 . M A E X , h = 1 m ∑ i = 1 m   h x i − y i , (21) R M S E X , h = 1 m ∑ i = 1 m   h x i − y i 2 , (22) MBE = 1 n ∑ i = 1 n   y i − y ^ i (23) R 2 = 1 − ∑ i = 1 n   h x i − y i 2 ∑ i = 1 n   h x ¯ i − y i 2 (24)

The Levenberg–Marquardt-BP neural network can accelerate the training process and is suitable for complex nonlinear regression problems. The goal of this algorithm is to adjust the network weight by minimizing the error function. The algorithm implementation process is as follows: 1) the gradient of the loss function is calculated. 2) The weights are updated according to the Levenberg–Marquardt algorithm. 3) Iterations are repeated until convergence or the maximum number of iterations is reached. The model calculation results are given below, where the forecast results are colored yellow, the drift data are colored blue, and the correct data are colored green. The performance of LM-BP-based IMU estimation is shown in Figure 8.

Gaussian process regression is a powerful non-parametric Bayesian regression method. It uses kernel functions to capture complex relationships of data by mapping the input space to a high-dimensional feature space. The cluster center of K-means (KM) is selected as the input feature point. The similarity between these feature points is calculated using the RBF kernel function. The model parameters are optimized by maximizing the posterior probability or minimizing the negative log-likelihood. The performance of KM-RBF-GPR-based IMU estimation is shown in Figure 9.

Two experiments are conducted, including the PSO–LSTM- and the PSO–ML-LSTM-based IMU estimations. The performances of both models are shown in Figures 10, 11.

The results of the experiments are shown in Table 3. Although KM-RBF-GPR slightly outperforms LM-BP in terms of RMSE and MAE, it exhibits a larger bias (MBE) and relatively lower goodness of fit. The performance of PSO–LSTM is superior to that of KM-RBF-GPR, particularly showing an improvement in the goodness of fit (R ²), but it still cannot surpass that of PSO–ML-LSTM. The RMSE of PSO–ML-LSTM is 8.3116, which is 1.9 units lower than that of LM-BP and 1.5 units lower than the RMSE of PSO–LSTM. The MAE and R ² of PSO–ML-LSTM are 6.1617 and 0.8594, indicating that PSO–ML-LSTM has no significant systematic bias and fits the data well. The average cumulative error of PSO–ML-LSTM is approximately 2.8°, which is significantly better than that of the other three models.

The ablation studies are conducted to evaluate the effects of the Kalman filter and random-forest feature selection. The results are shown in the Figures 12, 13. Compared with the result shown in figure, the results shown in Figure 8 indicate that the use of random-forest feature selection and Kalman filter can achieve a better performance of RMSE, MAE, MBE and R ².

A manipulator teleoperation experiment was conducted. By comparing the arm angles of the manipulator before and after IMU error compensation, the effectiveness of the proposed method was verified. The manipulator used was the 7-DOF manipulator described in the previous section, whose structure is shown in Figure 14. Figure 15 shows the arm angle conditions after 20 min, including 1) the manipulator’s arm angle without compensation (green curve), 2) the arm angle of the human operator (yellow curve), and 3) the manipulator’s arm angle using the proposed method (red curve). Without IMU error correction, the deviation between the manipulator’s arm angle and the operator’s arm angle is approximately 30°, and the deviation exhibits a nonlinear variation. After applying the proposed approach, the arm angle deviation is significantly reduced. The maximum deviation decreased from 30° to 10°, and the variation trends followed similar patterns.

Figure 16 shows the robot end-effector position, including 1) the actual end-effector position without compensation (green curve), 2) the operator wrist position (blue curve), and 3) the predicted position of the end-effector using the proposed approach (purple curve). Without IMU error correction, the position deviation between the end-effector position and the operator’s wrist is approximately 0.25 m, and the angle of the end-effector differs significantly from that of the operator’s wrist. After applying the proposed approach, the position deviation is reduced to 0.13 m, and the end-effector’s angle closely matches the operator’s wrist.