The rock drilling robot is mainly composed of the wingspan platform, left and right drilling manipulator, as shown in Figure 1. The wingspan platform C_w can be raised to enhance the work area. The drilling manipulator, which is mounted at the end of the wingspan, is the main working device of the rock drilling robot. The drilling manipulator has three moving joints and six rotating joints. The left and right manipulator oil cylinders C_d are placed at the front end of the manipulator to form a double triangle structure. At the same time, the left and right pitching cylinders C_s form a small triangle structure and are connected in series with the manipulator front support cylinder C_d to realize parallel linkage at the rear end of the manipulator. The extension of the manipulator telescopic cylinder C_b determines the longitudinal working distance of the manipulator. The hydraulic motor C_f of the flip joint is used to control the propeller around the axis Linear lateral rotation, with thruster tilt cylinder C_q to achieve different drilling angles. Compensation cylinder C_t pushes the bracket to move forward and backward to realize the fine adjustment of the position and posture. Each of the joints is actuated by hydraulic cylinders, which effectively serve as velocity sources.

DH method is a standard method to realize the transformation between joint variables and Cartesian coordinates and learn the kinematics modeling of the robot. However, when the joint axes between adjacent links are close to parallel, the position of the actual normal will deviate significantly from the theoretical normal due to the slight deviation of parallelism, which will affect the kinematic calculation.

Therefore, the MDH method is adopted to describe the conversion relationship of parallel joint axis systems. Its characteristic role is to add the rotation change Rot (y_i, β_i) around the y axis on the original DH model and set the initial value of β_i of adjacent joints to 0.

According to the structural characteristics of the rock drilling manipulator, the coordinate system transformation between links is designed, and the kinematic model of the manipulator is established, as shown in Figure 2.

When the coordinate systems of the manipulator 1, 2, 4, and 5 are parallel, the forward kinematics model of the rock drilling robot is established according to MDH method as shown in Table 1, and the homogeneous transformation matrix Tii-1 between adjacent link coordinate systems is obtained, as shown in Equation (1):

where m₁₁ = cθ_icβ_i, m₁₂ = −sθ_i, m₁₃ = cθ_isβ_i, m₁₄ = a_i−1, m₂₁ = sθ_icα_i−1cβ_i + sα_i−1sβ_i, m₂₂ = cθ_icα_i−1, m₂₃ = sθ_icα_i−1sβ_i−sα_i−1cβ_i, m₂₄ = −d_isα_i−1, m₃₁ = sθ_isα_i−1cβ_i−cα_i−1sβ_i, m₃₂ = cθ_isα_i−1, m₃₃ = sθ_isα_i−1sβ_i + cα_i−1cβ_i, m₃₄ = d_icα_i−1. s presents the sine function while c is the cosine function.

When the coordinate systems of two adjacent joints are not parallel, β_i = 0, and the second transformation Tii-1 between the links of two adjacent non-parallel joints is defined in Equation (2):

Substituting the parameters of each joint into the homogeneous transformation matrix, the transformation matrix of the robot manipulator end actuator relative to the base is obtained, as shown in Equation (3):

The geometric error mainly refers to the error between the ideal kinematic model parameters and the parameters of the actual MDH modeling (Khan and Chen, 2011; Cui et al., 2012). Therefore, the geometric errors are composed of measurement errors of length and angles during modeling, the calibration errors of the sensor, the manufacturing errors, and assembly errors. For positive kinematics parameters, a_i represents the distance between the joint axes of two adjacent joints, and the main error caused by a_i is the actual machining error. The symbol α_i represents the relative rotation angle of the two connected joints, and its main error comes from the coaxiality error in the assembly. d_i represents the relative distance between the two joints on the common joint axis, which is mainly derived from the measurement error during modeling. θ_i represents the rotation angle of two adjacent joints around the common axis; the error mainly comes from the observation error of the sensors.

The key role of the MDH method is to determine the axis of each joint. Therefore, by tracking and recording the trajectory of each joint's movement, the actual joint axis position of the robot can be established in the MDH model. The kinematic parameters of the manipulator are defined by error computation or estimation applying the posture measurement and position data of the end effector. In this article, the kinematic parameters are estimated by applying a differential evolution algorithm. The identified kinematics error parameters Δα, Δa, Δd, Δθ, Δβ are defined as: Kinematic parameters are estimated applying differential evolution algorithm. The identified kinematics error parameters Δα, Δa, Δd, Δθ, Δβ are defined as:

where e_d denotes the error of the end position, which indicates the distance deviation between the predicted position of the model and the actual position. The symbol e_δ is the error of the end attitude, which indicates the predicted rotation angle deviation between the model and the actual attitude on the x, y, and z axes. The observation prisms “r” and “s” are set at the front and end of the manipulator's compensation cylinder. The displacement from the deflection deformation and other nonlinear factors are ignored in the geometric compensation. Then the attitude error at the end of the robotic manipulator [edreδr] can be obtained from the position error of the two observation [edreds]. The geometric error model can be expressed as:

The differential evolution algorithm is a heuristic algorithm that can be applied to solve the optimization problem through competition and cross-mutation. In order to identify the parameters of the rock drilling robot, the common differential evolution algorithm is proposed as the following five steps:

Step 1 Initialization

The initial population is generated by randomly sampling the feasible search space defined by joint bounds, as shown in Equation (6):

Let n be the individual dimension, expressed as the value of the first individual of the generation population in the dimension, and the individual randomly generates an initial value that meets the constraints in each dimension:

In Equation (6), U_j and L_j represent the upper and lower bounds of the individual on the dimension and rand_ij(0, 1) represents a random number between [0, 1] that obeys a uniform distribution.

Step 2 Mutation

There are multiple solutions for the parameter identification of the rock drilling robot error model, the robust DE/rand/1/bin evolution model is used to preserve the diversity of the search population (Mlakar et al., 2015). The specific operation is to select an individual with the highest fitness from the current population, randomly select three different individuals, and scale the difference between the two vectors and add it to another individual vector to obtain a new variable Individual (Pavelski et al., 2016):

where F is the scaling factor of DE, the value range is [0, 1].

Step 3 Boundary constraint

The individual vector after the mutation operation may exceed the boundary of the search solution space, and the infeasible solution is converted to generate a new vector (Tian et al., 2020):

Step 4 Crossover

The differential evolution algorithm crosses the reference vector and mutation vector to increase the diversity of the next generation population. The specific operations are as follows:

Where: rand[0,1] represents the uniformly distributed random number on the jth dimension, C_r represents the crossover probability in the range of [0,1], j_rand represents a random number in 1,2,…,n.

Step 5 Selection

Based on the greedy selection mechanism, it compares the fitness of the mutated and cross-generated individual with the original individual and retains the highly adaptive individual as the next-generation population individual. The specific expression is as follows:

where f is the fitness function.

Since the rock drilling robot has three moving joints and six rotating joints, the dimension numbers of the error compensation model is 38. The high-dimensional solution requires a high number of populations and evolutionary algebra, which poses challenges to computer performance, algorithm model convergence, and computing speed. In order to accelerate the model convergence and increase the calculation speed under the same computing power input, a parallel structure of the differential evolution algorithm is proposed. The computing performance of the CPU can be enhanced by applying a parallel differential evolution algorithm to calculate independent individual populations. Specifically, due to the independence of the sub-populations, multiple processes can be allocated through different computing cores of the CPU to calculate the fitness of multiple sub-populations parallel. The structure of the parallel differential evolution algorithm is shown in Figures 2–6. Set the kinematic parameter error compensation vector Δα, Δa, Δd, Δθ, Δβ to be the population individual in the differential evolution algorithm. The sum of squared errors between the trajectory of the end of the robotic manipulator and the observation trajectory in the population is F.

Due to the heavy load, large span, redundant joints, and hydraulic character of the rock drilling robot, the error shows a high degree of non-linearity. After geometric error compensation is performed on the positive kinematics model, there are still errors due to the non-linearity. The non-geometric errors are composed of regular errors and irregular errors. Regular errors refer to predictable errors, including elastic deformation due to changes in load, frictional force direction, structural changes, and the cumulative error caused by the transducer. The irregular errors mainly come from observation errors and system errors. The observation errors are related to the dynamic position capture accuracy of the total station, the observation elevation angle, and the measurement accuracy of the optical sensor; the system errors mainly come from the command error of the encoder and the signal transmission communication delay, vibration errors caused by kinetic energy changes, etc. Regular non-geometric errors are related to the time and space of the robot arm movement. Take the compensation joint L3 as an example. Moreover, changes in the relative speed and direction of joint motion affect the magnitude and direction of frictional resistance, and thereby the position error is changed.

Features are important conditions that affects the machine learning models. Important features are extracted from the original data for the use of algorithms and models (Yu et al., 2017; Wen et al., 2021). Different algorithms have different feature selections during supervised learning training due to their training methods and structural characteristics. To perform supervised learning training on non-geometric error compensation of rock drilling robots, various features are extracted according to different input feature vectors. The feature types are divided into three types:

RBFNN is a three-layer neural network that can quickly approximate the nonlinear function to overcome local minimums. There are three layers in RBFN, namely including an input layer, a hidden layer, and an output layer, as shown in Figure 3.

The relationship of the RBFN from vector input to output is described as:

where x is the input vector, ω is the output layer connection weight, ϕ is the hidden layer basis function, and c_k is the basis function center of the kth neuron node

The radial basis kernel function is defined as Gaussian radial basis function:

where σ_k is the expansion width of the radial basis function of the first neuron in the hidden layer.

Features are composed of the joint sensor value and the three-dimensional coordinates calculated through positive kinematics. In this article, the input of the RBFN model is normalized “joint sensor value (9-dimensional)” and “positive kinematics model parameters (27-dimensional)” which is 36-dimensional in total, and the output is the predicted three-dimensional coordinates. When k = 100 and σ = 0.05, it is validated that the model shows the best result through experiments.

LightGBM is a fast, distributed, and high-performing framework based on gradient boosting decision tree (GBDT), which has wide applications in for decades for due to its superior predictive performance. Combining gradient boosting and decision tree, lightGBM has a strong expressive ability, good training effect, and can effectively prevent overfitting by controlling the growth of the tree. It can solve continuous and discrete values flexibly and is robust to outliers. The process is mainly composed of three steps, outlined as below:

Step 1 The training data set is given as: P = {(x₁, y₁), …, (x₂, y₂)}, let the loss function be L(y, γ), the learner is initialized:

Step 2 M base learners are iteratively generated

The negative gradient of each samples i = 1, 2, …, N is calculated as:

(x_i, r_mi)(i = 1, 2, …, N) is utilized as the training data for the next regression tree. The leaf node area R_(jm) of the m-th tree f_m(x), j = 1, 2, …, J,J is the number of leaf nodes of the regression tree. The best fit is calculated as

Step 3 The strong learner is generated as:

The endpoint error compensation is used as the training label for supervised learning; its training features have a total of 66 dimensions, which are “feature of the joint (9 dimensions)” + “feature of positive kinematics model (27 dimensions)” + “feature of control system (30 dimensions).” K-fold cross-validation method is applied. The test data D is divided into k parts in equal proportions. One of the parts is applied as the verification set. The other K−1 parts of data are used as the training data. After performing K times of training, the average result of the K training experiments is used as the final result. When the model stability is low, increasing the value of K can achieve better results.

In order to achieve accurate position compensation for the rock drilling robot, the error compensation experiment is performed. In this experiment, the SWDT82 rock drilling robot developed by Sunward Intelligent Co., Ltd, as shown in Figure 5, has a maximum manipulator length of 12.8 m and a wingspan of 5.2 m. The maximum coverage section width of the dual-manipulator rock drilling robot is 17.65 m, the maximum coverage section height is 13.38 m, and the maximum invert depth is 3.31 m.

The rock drilling robot is equipped with a three length transducer and seven angle sensors as shown in Figure 6. The angle sensors are the rotary encoder of POSITAL's model IXARC, and its theoretical measurement error is 0.1°; the length transducer is TRANSTRONIC's model 022441, and the theoretical error of measurement linearity of this sensor is 0.1% mm. The observation instrument in the experiment is a Topcon laser total station, which has functions such as motor tracking, automatic sighting, and Bluetooth transmission.

To compensate for the geometric error of the rock drilling robot, the motion of a single joint is carried out. The position and posture of the rock drilling robot are collected through Topcon total station transmitting to the master computer using Bluetooth as shown in Figure 7. The data of the angle sensors and length transducers are transmitted to the master computer w, b, r, and s corresponding to the joint 3, the joint 5, the front of the compensation cylinder 9, and the end position of the compensation cylinder 10, as shown in Figure 8. Moreover, the end posture of the rock drilling manipulator is determined by the position of observation point r and s.

So as to identify the geometric error of the manipulator, single-joint motion of the manipulator is carried out. After substituting each observation point, obtained from sensors, into the forward kinematics model for calculation, the absolute errors can be obtained. The absolute average error of each observation point is er0̄: es0̄: eb0̄: ew0̄ = 323.59: 477.53: 895.52: 930.88 (mm).

Afterwards, the multi-joint motion experiment is carried out. Trajectory data of the observation points of rock drilling robot motion r, s, b, w is applied for error analysis. In the multi-joint movement, the absolute average error of each observation point data set er0̄: es0̄: eb0̄: ew0̄ = 355.61: 783.17: 1085.32: 1288.916 (mm).

In order to identify the forward kinematics of the manipulator, the differential evolution method is applied to compensate for the geometric error. This experiment not only explores the error compensation effect of the error model method but also explores the influence of the amount of data on the convergence of the parameter identification of the error model method. Due to the high-dimension of the model, the population of the differential evolution algorithm is set to 30,000. To avoid premature convergence of parameter optimization, the differential cross recombination probability C_r is set to 1, the scaling factor F = 0.5, and the number of evolutions is set to 400. To improve the performance of the model, the data set is shuffled and substituted into the model multiple times to obtain the optimal solution. To optimize the error model, the joint constraints are increased. Set the end trajectory of each observation point as Lw′, Lb′, Ls′, Lr′, the trajectory observation point set is P′ and the model prediction point set are (X′, Y′, Z′, Then the position error of each observation position are ewd,ebd,esd,erd. The average absolute error F of the optimal error model with ideal multi-joint constraints at each observation point is expressed as:

The solutions of F_w, F_b, F_s, and F_r are difficult to satisfy simultaneously, and the reduction of the error of some joint clusters requires the degradation of the error fitting effect of the other joint clusters, resulting in the Pareto solution. And new observation errors are introduced by adding trajectory constraints, which reduce the overall accuracy of the model. In order to solve the above problems, an equal amount of data in the model calculation are applied for each joint. Moreover, a weight coefficient is assigned to each joint, and the functions are obtained by linearly combining the multi-objective functions:

where μ_r, μ_s, μ_b, and μ_w are the weight coefficients of the multi-objective error constraint, andn_r, n_s, n_b, and n_w are the numbers of each data set. The parameters of forward kinematics of the optimal method are obtained shown in as Table 2.

The average absolute error can be calculated from the compensated positive kinematics. A parallel differential evolution algorithm is applied to move a single joint as the first experiment. Then, a parallel differential evolution algorithm is applied for moving multi-joints applying as the second experiment. Ultimately, the optimal differential evolution algorithm is applied as the third experiment. Consequently, the distribution of errors is shown in Table 3.

Applying the optimal method of parallel differential evolution, the errors of the end effectors are compensated, and the best performance is achieved compared to the parallel differential evolution algorithm.

The overall control flowchart is shown in Figure 9. Firstly, the modified positive kinematics is obtained through an optimal parallel differential evolution algorithm to compensate for the geometric error. Afterward, the features are applied to train RBFN. The estimated position and features are combined to train light GBM. The non-geometric errors are compensated through RBFN and light GBM. The top 10 features are shown in Figure 10. X_i, Y_i, Z_i (i=1, 2..) are the Cartesian coordinates of the joint i, orient3 is the control current of the hydraulic valve, Propeller_Flip_Angle is the sensor of the flip joint, Z_dabi, Y_dabi denote the Z coordinate and Y coordinate of the big manipulator separately, Y coordinate. The hydraulic valve current control ranks 7th in feature importance as a system feature, indicating that the system feature has an important contribution for the model.

The absolute error of the endpoint is obtained for multiple experiments in Figure 11. All of the experiments are carried out by applying the same positive kinematics from the optimal parallel differential evolution. The average error of the RBFN is 41.5 mm, while the standard error of the RBFN is 33.5 mm. Moreover, the average error of the lightGBM is 41.0 mm, while the standard error of the RBFN is 28.5 mm. Moreover, the average error of the hybrid method is 30.5 mm, while the standard error of the hybrid method is 17.4 mm. It can be concluded that the hybrid method can achieve a good performance.

A validation experiment is carried out to verify the effectiveness of the proposed method as shown in Figure 12. Two hundred points collected along the trajectory are analyzed in two error compensation methods. The blue line is the absolute error before compensation. The yellow line is the absolute error applying the parallel difference algorithm for geometric compensation. The green line is the absolute error applying RBFN combined with light GBM for geometric compensation and non-geometric compensation. It can be concluded that the error of the rock drilling robot is compensated effectively by applying the proposed method, which meets the practical application needs of rock drilling robots.