Kinematic Calibration of Industrial Robots Based on Distance Information Using a Hybrid Identification Method

To improve the positioning accuracy of industrial robots and avoid using the coordinates of the end effector, a novel kinematic calibration method based on the distance information is proposed. (e kinematic model of an industrial robot is established. (e relationship between the moving distance of the end effector and the kinematic parameters is analyzed. Based on the results of the analysis and the kinematic model of the robot, the error model with displacements as the reference is built, which is linearized for the convenience of the following identification. (e singular value decomposition (SVD) is used to eliminate the redundant parameters of the error model. To solve the problem that traditional optimization algorithms are easily affected by data noise in high dimension identification, a novel extended Kalman filter (EKF) and regularized particle filter (RPF) hybrid identification method is presented. EKF is used in the preidentification of the linearized error model. With the preidentification results as the initial parameters, RPF is used to identify the kinematic parameters of the linearized error model. Simulations are carried out to validate the effectiveness of the proposed method, which shows that the method can identify the error of the parameters and after compensation the accuracy of the robot is improved.


Introduction
With high generality and industrial flexibility, robots have been widely used in modern manufacturing, which play an important role in the fields of automobile manufacturing, logistics, machinery manufacture, and so forth. Traditionally, industrial robots are used in the way of teaching programming, for example, in spot welding, handling, and other similar repetitive works. In recent years, the position accuracy of industrial robots is required to be higher and higher in the fields of mechanical manufacturing, carving, measuring, testing, and so forth.
About 90% position errors are caused by the inaccuracy of the kinematic parameters in the controller [1]. e main way to improve the position accuracy of industrial robots is kinematic calibration by which the structural parameters errors due to the tolerance in machining and assembly of robots can be identified and compensated [2]. ere are two levels for kinematic calibration. e first level is to identify and compensate the errors between the transducer reading of joints and the actual joint angle. e second level is to identify and compensate all the kinematic parameters of robots. Generally, the kinematic calibration of robots can be classified into four steps [3]: modelling, measurement, parameter identification, and compensation. e kinematic parameter identification of industrial robots is a high nonlinear problem. To obtain precision identified parameters, the kinematic model is required with the features of completeness, continuity, and minimality [4]. e standard Denavit-Hartenberg (DH) model [5] is the most commonly used model in kinematic modelling of serial open-chain robots. However, there will be singularity when two adjacent joints are completely or almost parallel. Hayati [6] proposed a modified DH (MDH) method, which introduces a rotation parameter to solve the singularity problem. S-model was also proposed [7] to solve the singularity problem by Stone. A Complete and Parametrically Continuous (CPC) kinematic model was proposed by [8] to provide a kinematic modelling method with completeness and continuity.
In the step of measurement, a group of poses and positions of the robot are acquired for the following parameter identification usually by laser trackers [9][10][11][12] or other coordinate measuring instruments [13][14][15]. Nadia Schillre [16] reported a calibration method with a laser tracker to acquire data. e maximum position error of the robot is reduced from 0.4883 mm to 0.1801 mm. Albert Nubiola [17] presented a calibration case for ABB IRB 1600 using laser tracker, in which the mean position error of the robot was reduced from 0.968 mm to 0.696 mm. With the advantages of high precision, dynamic tracking, and noncontact measuring, laser trackers become the main measurement instrument in the calibration field of robots. But they are so expensive that ordinary users cannot afford them. Except laser trackers, coordinate measuring machines (CMMs) [18], measuring arms [19], ballbars [20], and so forth are also used in calibrating robots. Actually, the coordinates of the end effector of the robots are acquired by instruments in the previous researches [21][22][23], which are needed in the following identification. ese coordinates are acquired in the measurement coordinate system of the instruments, which have to be transformed into the base coordinate system of the robot, and this transformation will bring error and increase the difficulty of identification. Unlike coordinates, distance is a relative value, which is always the same in different coordinate systems without any transformation. e third step, identification, actually is an optimization problem, in which the least square (LS) algorithm [24][25][26][27] is the most commonly used way. But LS is easily affected by the noise of data, which can lead to failure in identification. Extended Kalman filter (EKF) algorithm was more suitable for the identification of the robot, the working principle of which is similar to that of Kalman filter [28]. e expected EKF deals with nonlinear systems that retain the Taylor expansion first-order terms of the nonlinear function and ignore the higher-order terms. But EKF cannot deal with high nonlinear system and only suit Gaussian noise system [29]. Actually, the kinematic model of the robot is a highly nonlinear system.
Referring to ISO-9258 [30,31], the distance accuracy of industrial robot can be defined as the deviation between the measured distance and the command distance between the continuous moving points in workspace. To avoid using coordinate measuring instruments and the transformation between the measurement coordinate system and the base coordinate system of the robot, we propose a kinematic calibration method based on distance error modelling and the definition of distance accuracy to improve the position accuracy of the robot. A regularized Particle Filter (RPF) [32] is used to improve the ability of identification for the kinematic model of the robot in the light of the excellent performance of sample filtering method in solving highly nonlinear problem in identification. Since RPF is sensitive to its initial value, EKF is used in preidentification of the kinematic parameters of the robot. e results of the preidentification are used as the initial value of RPF. e paper is organized as follows: In Section 2, we introduce the kinematic modelling and distance error modelling. In Section 3, the identification of the robot is discussed and the redundant parameters in the identification model are analyzed by numerical method based on distance error. Comparative simulations are conducted in Section 4, and some conclusions are given in Section 5.

Kinematic Modelling and Distance
Error Modelling

Kinematic
Modelling. e controller of industrial robots is essentially a semi-closed-loop system, which sends control instructions to the joint servo motors and acquires the positions of them. e positions of the joints are calculated from the positions of the servo motors, which are used to determine the poses of the industrial robots by kinematic model.
An MDH method [33] is used to build the kinematic model of the industrial robot, in which the relationship between two adjacent rotation joints is described by a homogeneous transformation matrix. e kinematic model consists of the controller, which consists of four groups of parameters: link length a, link offset d, link twist angle α, and joint angle θ [33]. First, the coordinate systems of a 6-degreeof-freedom (DOF) industrial robot named ER20-C10 are established according to the MDH method, as shown in Figure 1. en, the nominal values of kinematic parameters of the robot are obtained, as shown in Table 1.
According to MDH method, the transformation matrix T i,i−1 of two arbitrary adjacent rotation links {i−1} and {i} is as follows: For the generic industrial robot with N t -DOF, the transformation matrixes from the base frame {0} to the end effector can be written as where T K,E � 0 0 l z 1 T and l z is the end effector (EE) offset along the axis of joint N t . e first three rows of the last column of T 0,E denote the position of EE, that is, P(·).

Distance Error Modelling.
According to ISO-9258 [30,31], between the continuous moving points during the movement of robot in workspace, the distance accuracy of industrial robot refers to the deviation between the measured distance and the command distance, illustrated as Figure 2. e distance of the two adjacent points (i and i + 1) in workspace can be expanded as e position error of the ith point of the acquisition data is defined as where J is the Jacobean matrix and X is the error of kinematic parameters. According to equations (3) and (4), one can obtain the distance error formula shown as follows: Equation (5) indicates that the distance error can be expressed by partial pose information. Furthermore, equation (5) can be rewritten as [34] Δl · P e (i + 1) − P e (i).

(6)
According to equation (6), the measured distance can be calculated by two corresponding adjacent points [34]. erefore, at every acquisition position, one has to measure twice to obtain the distance used, which is time-consuming and labor-consuming. In this paper, we proposed a modified distance model, which is more convenient and less timeconsuming than the general one. e original position of the device is set as P 0 (x p0 ,y p0 ,z p0 ) in the frame of the measurement device, which is fixed. P 0 is considered to be error-free, which is also set as one of the two adjacent points in equation (6). Hence, equation (6) can be rewritten as where l m (i) is the ith distance from EE to P 0 , which can be directly measured by the measurement system. P 0 is known as the fixed point of the measurement system. e original error of the measurement system can be neglectable for the distance is a relative value. With P 0 being the fixed reference point, to get the distance, one only needs to measure one time at a position. us, the measuring process can be simplified with proposed strategy. Equation (7) describes the relationship between the distance error model and the  No. of joints Complexity position error. Assuming that the deviation of the parameters is small, (7) can be linearized as the first-order term of its Taylor expansion around the nominal values of the kinematic parameters, as shown in the following equation: where J ∈ R 3n×N , n is the number of calibration distances, and N is the number of parameters. e Jacobian matrix J in distance error model is To define the problem of identification with filtering algorithm, the errors for kinematic parameters X and the measured data Y are described in the state space, which can be rewritten as follows: X � Δα 1 · · · Δα 6 Δa 1 · · · Δa 6 Δθ 1 · · · Δθ 6 Δd 1 · · · Δd 6 T , where k is the number of iterations and X k represents the parameters' errors. Y k is the corresponding calibration distance measured by the measurement device, u k is the process noise of measurement, and h (·) is the measurement function.

Kinematic Parameter Identification
e effect of gradual error accumulation of links decreases the absolute position accuracy significantly for the openchain structure of industrial robots. Kinematic parameter identification is an important way to improve the absolute position accuracy.
In the kinematic parameter identification, the data of the joint angles and the corresponding distances are required. e forward kinematic solutions of the robot can be obtained by combining a group of joint angles, with which the distances from EE of the robot to P 0 can be calculated. Unlike the kinematic parameter identification using laser tracker to acquire the coordinates of EE, the proposed method in this paper only needs the distances acquired by pull wire sensors or another distance measurement system.

Preidentification Based on EKF.
e solution for function (8) can be obtained commonly using a nonlinear LS approach. Although nonlinear LS approach is a fast and efficient algorithm that can be used to solve nonlinear equations, it is very sensitive to noise. To obtain more accurate and stable solutions, EKF algorithm is used in preliminary identification.
e EKF algorithm consists of two steps: prediction and update [29]. e parameters' error X can be regarded as the state vector to be preliminarily identified. First, the prior estimate of state vector X k+1 and the covariance matrix are predicted by the former state vector X k and P k .
where k denotes the number of iterations, X k is N-state vector consisting of parameters' deviation at kth iteration, X k+1 is the prior estimation of state vector, P k is the covariance matrix of the estimation state, and Q k is an N×N covariance matrix of the system process noise.
With the prior estimate of state vector, the optimal Kalman gain can be calculated with (12), and the covariance matrix of state vector is obtained, as shown in (14). en, the state vector can be updated by the measurement data with (13).
where K k + 1 is the Kalman gain at kth iteration and R k + 1 denotes the covariance matrix of measurement noise. Y k + 1 is the distance error measured by the measurement system. J is the Jacobian matrix. X k + 1 is the posterior estimate of state vector. e details of the EKF algorithm are illustrated in Figure 3.

e RPF Algorithm.
When the noise of the acquired data is non-Gaussian, it is difficult for the EKF algorithm to get a reliable identified result. e RPF algorithm does not need the model to be linear or compellent hypothesis for Gaussian noise [35], the basis of which is particle filter. e preidentification result of EKF is applied as the initial value of the RPF which can identify the state of nonlinear and non-Gaussian noise system. With a certain amount of random state samples (i.e., particles), the RPF algorithm can approximate the posterior probability density function of the identifying model. It provides a suboptimal solution with a finite number of samples. All state samples are transferred in dynamic systems through importance sampling. Subsequently, it updates the posterior distribution sequentially.
Apart from the resampling step, the RPF is similar to generic sampling importance resampling (SIR) filter [36]. Considering the robot that is described by the state space model, the RPF algorithm operates N p particles to approximate the posterior density. e RPF algorithm also consists of two steps, prediction and update, similar to EKF. But it is described statistically with the state transition probability density p(x k |x k−1 ) and the observation likelihood probability density of the state vector p(x k |y 1:k ).
In the prediction step, the predicted state vector is estimated by the last time state X k-1 and its measurement data With the latest measurement data Y 1:k , the predicted state can be updated; then, its approximated discrete posterior density can be obtained: where p(Y k |Y 1:k−1 ) is defined as the normalization constant, as shown in (20).
With Monte Carlo method [37], the approximated discrete posterior density can be simplified, as shown in the following equation: where δ(·) is the Dirac delta function, N p is the number of particles, and ω i k is the weight of the ith particle at kth iteration. p(X k |Y 1:k ) is the approximated discrete posterior density. According to Bayesian estimation [37], all the particles need to sample from the target posterior density with the knowledge of observation Y 1:k , but this is unrealistic. Importance sampling [37] is introduced in the RPF algorithm, as well as the importance samples from a known posterior density (i.e., importance density). In the identification, the weight of the particles and the normalized weight can be defined as (19) and (20), respectively.
where e d k is the distance error of the ith particle. According to Bayesian estimation [37], the minimum mean square (MMS) estimate X i k can be described as the expectation of identified state vector combining all particles [38], as shown in the following equation: In the identification model with distance error, the weight of each particle is evaluated using its distance error. e particles with lower distance error Δl have higher weight; that is, After a number of iterations, the weight will occupy on one particle (almost close to one) and the weights of other particles are so small that they will be negligible due to the degeneracy of particles [37]. e degeneracy problem indicates that abundant computation is applied to updating particles whose weight is almost zero. e negligible particles with tiny weight make no contribution to approximation of the target posterior density. Resampling [39] should be applied to avoid this problem. Note that an effective sample size N eff is introduced to measure the degeneracy of particles [39]. If N eff is too small, resample will be triggered.  Complexity Resampling is an approach to prevent particles from degeneracy, which will introduce the decay of diversity of particles in turn. is problem leads to poor approximation of the posterior density. One way to solve this problem is modifying the resampling, in which the particles are resampled from the continuous approximation of the posterior density [40]. e approximated continuous posterior density in the RPF [40] is where where k h (·) is the Kernel probability density function, h > 0 is the Kernel bandwidth, and n is the dimension of the estimated state vector X. k h will be used for resampling particles on continuous approximated density. With appropriate function k h and h, the mean integrated square error (MISE) between real and regularized empirical posterior densities in (24) can be minimized [38]. e details of the RPF algorithm are illustrated in Figure 4. It is assumed that the coordinate of P 0 is constant, which can be regarded as one of identified parameters. e offset of the adapter of the measurement device should also be considered. In general, a complete kinematic model includes 28 unknown kinematic parameters, which need to be identified in our framework. e distances used in identification should be more than 10 mm to increase the certainty of identified results. Figure 5 shows the flow chart of the identification.

Determination of Redundant Parameters.
Singularity analysis of the identified matrix J is very important in identification calculation, which can lead to nonconvergence or inaccurate results. Generally speaking, the maximum number of parameters is 4r 1 + 2p 1 + 6 in a generic serial-link robot with r 1 rotation axes and p 1 translation axes [41]. To obtain reliable results with the identification algorithm, redundant parameters must be determined and excluded from the error model before identification. e condition number of the identified matrix is where κ(·) is the measure of sensitivity or stability of the identified matrix with regard to a small change in the input argument. In high dimension full-parameter identification, with a larger condition number, the identified matrix J is often nonfull rank. In this paper, the singular value decomposition (SVD) [42] is applied to reduce the condition number and avoid singularity problem. Considering total Jacobian matrix H � [J 1 , J 2 , . . ., J n ], it can be written as where U ∈ R 3n×3n , V ∈ R 25×25 , and they are orthogonal matrices; H ∈ R 3n×N , Σ � diag(σ 1 , σ 2 , σ 3 , . . . , σ r ), r is the rank of H, n is the number of calibration distances, N is the number of parameters, and, in this case, N is 25. Hence, the number of the parameters that are redundant can be obtained as 25 − r.
After SVD calculation, r can be determined, that is, r � 20. erefore, there are 5 redundant parameters that need to be eliminated from the process of identification. According to the property of distance error, the original coordinates of P 0 are determined in the base frame. From Figure 1, it can be found that the four parameter errors in joint 1 may cause the overall change of the global coordinate system. Hence, P 0 will change accordingly, which makes it identifiable. erefore, there are 19 parameters including P 0 (x p0 , y p0 , z p0 ) which need to be identified; and they can be written in the format of a vector: X � Δα 1 · · · Δα 5 Δa 1 · · · Δa 4 Δθ 2 · · · Δθ 5 Δd 2 Δd 4 Δd 6 T .

Verification
To verify the proposed kinematic parameter identification method, a 6-DOF industrial robot ER20-C10 was used, whose maximum load and repeatability are 20 kg and 0.08 mm, respectively. e nominal values of parameters are shown in Table 1. e settings of the error for each parameter are shown in Table 2. It should be noted that the error of the redundant parameters is set as zero, since these parameters are not independent. e two termination conditions of the iteration are as follows: the maximum number of iterations is 10000 and the minimum objective function is less than 10 −6 . Both of the covariance matrices Q and P were initialized as 10 −4 × I 3 * 3 in EKF; and the noise of the covariance matrices of measurement R is set as 10 −4 × I3 * 3 in the simulation. e bandwidth h is set as 0.42. A random and bounded measurement disturbance [−0.02 mm, 0.02 mm] is added to each distance to simulate the measuring error in the actual device.
To verify the effectiveness of the proposed algorithm, a group of 100 poses were generated randomly for identification and verification. Considering the workspace of the robot, the six joint angles were chosen randomly from To compare the effectiveness, three different algorithms (i.e., LS, EKF, and EKF + RPF) were used to identify the kinematic parameters that were compensated to the control model of the robot. e maximum distance error was used as the index to evaluate the results of compensation. Figure 6 is the distance error distribution chart of the robot after compensation by LS, which shows that all the distance errors are reduced greatly, and the average error is from 0.4827 mm to 0.1120 mm by a reduction rate of 76%. Figure 7 is the distance error distribution chart of the robot after compensation by EKF, which shows that the average error is reduced from 0.4827 mm to 0.1352 mm by a reduction rate Calculate ω k and total weight ω k i i iω k = (ω k ) -1 ω k    Complexity of 72%. Figure 8 is the distance error distribution chart of the robot after compensation by EKF + RPF, which shows that the average error is reduced from 0.4827 mm to 0.0780 mm by a reduction rate of 84%. Table 3 shows the statistics of the maximum, mean, and RMS values of the distance errors before and after compensation by the three algorithms. By comparison, one can find that the effectiveness of RPF is the best of the three.

Conclusion
is paper proposed a novel kinematic calibration method based on the distance information for the industrial robot, in which an EKF and RPF hybrid identification method was adopted. To simplify the data acquisition procedure and also avoid using expensive coordinate measuring instruments, the moving displacements of the end effector were used as the reference to build the kinematic parameter error identification model. en, the error model was linearized by differentiating of the coefficient matrix. By SVD of the coefficient matrix, the redundant parameters were determined and removed from the list of parameters to be identified. e preidentification of the parameters was conducted with EKF, by which we obtained a group of initial parameters of the kinematic model. RPF was used to identify the kinematic parameters of the linearized error model with the initial parameters. Comparative simulations were conducted based on the linearized kinematic error model. e compensation results of the three algorithms are all effective in improving the accuracy of the industrial robot, but the hybrid algorithm proposed in this paper is more precise and fast in iteration calculation. In the future work, precision distance measuring device will be developed, and effectiveness of the device and the proposed algorithm will be further verified by experimental tests.
Data Availability e data were curated by the authors and are available upon request.

Conflicts of Interest
e authors declare that there are no conflicts of interest regarding the publication of this paper.