# End-Effector Pose Estimation in Complex Environments Using Complementary Enhancement and Adaptive Fusion of Multisensor

Authors: Mingrui Luo; En Li; Rui Guo; Jiaxin Liu; Zize Liang
Journal: Journal of Sensors
 (2021)

---

## Abstract

Redundant manipulators are suitable for working in narrow and complex environments due to their flexibility. However, a large number of joints and long slender links make it hard to obtain the accurate end-effector pose of the redundant manipulator directly through the encoders. In this paper, a pose estimation method is proposed with the fusion of vision sensors, inertial sensors, and encoders. Firstly, according to the complementary characteristics of each measurement unit in the sensors, the original data is corrected and enhanced. Furthermore, an improved Kalman filter (KF) algorithm is adopted for data fusion by establishing the nonlinear motion prediction of the end-effector and the synchronization update model of the multirate sensors. Finally, the radial basis function (RBF) neural network is used to adaptively adjust the fusion parameters. It is verified in experiments that the proposed method achieves better performances on estimation error and update frequency than the original extended Kalman filter (EKF) and unscented Kalman filter (UKF) algorithm, especially in complex environments.


---

## Body

## 1. Introduction

In special manufacturing, troubleshooting, and other fields, there will inevitably be a closed and narrow operating environment, which is high risk and inefficient for manual operation. Therefore, the research on robots that can be used in this complex environment is significant and challenging, especially since the traditional manipulators are unable to do the job due to their bulky structure.Kinematic redundancy is one of the most important properties that determine whether the robots can accomplish the task in the highly constrained environments mentioned above [1]. Manipulators with kinematic redundancy are called redundant manipulators. Redundant manipulators can be implemented in different ways, such as the wheeled mobile manipulator used for part handling [2, 3], the hyperredundant manipulator used for engine manufacturing [4] or spacecraft maintenance [5], and the flexible bionic manipulator inspired by the octopus claw [6] or the elephant trunk [7]. Although redundant manipulators have excellent kinematic performance, they also face the challenge of control accuracy due to the accuracy of end-effector pose estimation.Traditional manipulators generally use joint encoders as sensor inputs and use kinematics models to calculate the pose of the end-effector as information feedback for motion control. Since the accuracy of the joint encoder can be very high, the pose calculated by the kinematics model has low noise. However, redundant manipulators usually have a large number of joints and long slender links, which means the pose will be affected by the transmission gap, mechanical vibration, elastic deformation, and other factors that deviate from the calculated value. As a result, the manipulators are in an uncertain state of kinematics and dynamics, so complex control models need to be established to ensure control accuracy [8–10].Complex control models may introduce system uncertainty, so many researchers try to solve the problem of unreliable encoders in redundant manipulators by studying sensor processing methods. Analyzing the source of the measurement error and establishing a compensation model can effectively improve the sensor accuracy. For example, in a cable-driven redundant manipulator, the accuracy of the encoder-based pose estimation can be improved by correcting the cable hole location error, cable length error [11], and nonnegligible cable mass [12].In addition to error correction for a single type of sensor, multisensor fusion is also the main research direction. The most commonly used sensors fused with encoders including angle sensors, vision sensors, and inertial sensors. For one thing, adding angle sensors at joints and fusing with motor encoders can correct transmission errors and improve pose estimation accuracy [13]. For another, adding one or more eye-to-hand cameras to the environment and performing the visual measurement on redundant manipulators can correct deformation errors [14] and is also an effective method for flexible robot pose estimation [15]. Moreover, the Kalman filter and its improved methods are also used to fuse motor encoders and inertial sensors, which have good dynamic response performance and can correct vibration errors [16]. Because the sensor and the fusion algorithm are lightweight, it is especially suitable for small bionic robots [17].The sensor processing methods mentioned above can indeed solve the problem of end-effector pose estimation, but there are also many unsatisfactory shortcomings, especially when applied in complex environments. Therefore, the estimator constructed in this paper uses the compensation model and multisensor fusion at the same time, but in different processing stages. In a dynamic environment, although the preestablished compensation model is prone to failure and cannot accurately reflect the relationship between the sensor and the estimated state, it can effectively correct the wrong measurement [18]. Therefore, this paper does not directly use the compensation model for state estimation but as the preprocessing part of the state estimation to enhance the original data of the sensor. In terms of multisensor fusion, although visual-inertial fusion has become an effective method of mobile robot navigation [19], it is not common in end-effector pose estimation of redundant robots, and it is mostly the fusion of eye-to-hand camera and kinematic sensors [20]. In a narrow environment, it is difficult to apply an eye-to-hand camera, so the vision sensor needs to be arranged at the end of the manipulator, which means that the vision sensor is more susceptible to the influence of obstacles. Therefore, based on the Kalman filter algorithm, this paper establishes a nonlinear motion prediction of the end-effector and a synchronization update model of the multirate sensors that combines vision sensors, inertial sensors, and encoders. Moreover, since the neural network can improve the robustness of pose estimation [21], this paper introduces the radial basis function neural network to adjust parameters adaptively. Through the above-mentioned improved fusion strategy, the pose estimation of the end-effector can be more adaptable to complex environments and can better overcome the adverse effects of sensor noise and different sampling rates of multisensor fusion.The rest of this paper is organized as follows: Section2 introduces the overall method. Section 3 introduces the complementary enhancement of the original data. Section 4 introduces the fusion strategy including the nonlinear prediction, the synchronization update, and the adaptive adjustment. Section 5 introduces the experimental results and analysis, and Section 6 summarizes the paper.


## 2. The Overall Method of the Pose Estimation

Figure1 shows the small redundant robot studied in this paper, which can enter the narrow space to accomplish the specified tasks.Figure 1
The small redundant robot (a) end-effector, (b) redundant manipulator, (c) base, and (d) connecting rod.Due to the miniaturized mechanical structure, when the end-effector is equipped with a load, the connecting rod between the joints will inevitably be bent and deformed. This deformation error will accumulate at the end-effector; as a result, the pose estimation based only on the joint encoders produces deviation. Therefore, structured-light RGB-D camera, magnetic angular rate, and gravity (MARG) sensor are added, together with the joint encoders.In this paper, the fusion pose estimation method proposed mainly includes two parts: preenhancement and adaptive fusion. The system block diagram is shown in Figure2.Figure 2
Block diagram of the pose estimation method.The RGB-D camera and the MARG sensor can compensate for the error caused by a single encoder measurement to a certain extent, but at the same time, it is easy to reduce the robustness of the system due to the interference of the external environment. Considering that the RGB-D camera and the MARG sensor themselves are integrated sensors composed of multiple sensing units, each sensing unit has complementary characteristics, which can be used to enhance the original data of the sensors so that the fusion of multiple sensors can be completed based on robust data. According to the above theory, preenhancement is added before multisensor fusion.The preenhancement part includes the image enhancement of the RGB-D camera and the quaternion extended Kalman filter (QEKF) fusion attitude calculation of the MARG sensor. In the image enhancement part, the color image is used to repair the depth image, which can make the measurement of the RGB-D camera more accurate and reliable in a narrow environment (Section3.1). In the QEKF part, the magnetometer, accelerometer, and gyroscope, which are all parts of the MARG sensor, are fused to correct the drift error of each independent sensing unit (Section 3.2).Through the preenhancement part, the robustness of each independent sensor is improved, but the following problems still need to be overcome in the multisensor fusion process: (a) nonlinear state model, (b) the sampling rate of different sensors differs greatly, and (c) sensor abnormal handling in extreme environments.In the adaptive fusion part, the RBF and adaptive UKF (RBF-AUKF) method is modified to solve the above problems, which mainly includes pose prediction, synchronization update, and adaptive adjustment. In the pose prediction part, by establishing the pose prediction model of the end-effector and the nonlinear processing method, the pose statexk∣k−1 at the next moment can be predicted based on the pose state xk−1 at the previous moment (Section 4.1). In the synchronization update part, a synchronization check is performed on the preenhanced quaternion attitudeqV and position pV obtained by the visual odometer; quaternion attitudeqM, angular velocity ωM, and acceleration p¨M obtained by the QEKF algorithm; and quaternion attitudeqE, position pE, and velocity ṗE obtained by the kinematic calculation. According to the check result, the control variable is relatively updated when the sensor data zk is asynchronous, while the state variable is absolutely updated when the sensor data zk is synchronous. After the prediction and update, the fusion pose estimation xk at this moment is obtained (Section 4.2). In the adaptive adjustment part, this paper selects the prediction error xk−xk∣k−1, estimation error covariance matrix Pk, Kalman gain Gk, and observation error zk−zk∣k−1 after updating as features. Based on the RBF neural network, the process, and observation noise covariance matrix Qk, Rk is adaptively adjusted to improve the robustness of fusion (Section 4.3).


## 3. Sensor Enhancement Based on Complementary Characteristics

### 3.1. Enhancement of the Depth Image

The structured-light RGB-D camera is composed of a color camera, an infrared camera, and an infrared transmitter and obtains the depth value through the speckle feature of the infrared structured light projected on the surface of the object [22]. In a closed and narrow environment, infrared light is easily reflected and absorbed by objects multiple times, so it is easier to produce invalid depth values locally, which in turn introduces additional noise for pose estimation.Since the invalid depth value occurs locally, the invalid point can be repaired by a valid point that is in a similar region to the invalid point. Since the color camera is not affected by the narrow environment, it can be used to determine whether the valid point and the invalid point belong to the same type of patch gathered by the point cloud, that is, the similar region. Therefore, the enhanced method proposed in this paper includes two steps, namely, similar region extraction and local point cloud reconstruction, as shown in Figure3. The former extracts similar regions in the depth image through the texture features of the registered color image. The latter uses valid points in the similar region to construct a local model of the region and then recalculates invalid points in the region with the model to complete the enhancement repair process.Figure 3
Depth image enhancement in a narrow environment.In the similar region extraction step, the Mean Shift algorithm is used to extract the color texture information of the color image, and the Flood Fill algorithm is used to extract the corresponding color-connected domains to form a similar region mask. Finally, the mask is applied to the depth image to obtain the region segmentation of the similar point cloud. Since the Mean Shift algorithm needs to perform multiple iterative calculations for each pixel, it will affect the real-time performance of the system in practical applications. In order to speed up the calculation process, this paper downsamples the input color image to a lower resolution, then forms the similar region mask on the low-resolution image, and finally performs bilinear interpolation on the mask to restore the original scale.In the local point cloud reconstruction step, by fitting the patch model in the similar region and recalculating the depth value of the invalid point on the patch, the invalid point is repaired. The least-squares method is the basic method to fit the patch model, but in practical applications, the region segmented by the color image has deviation, and not all valid points can be used to fit the patch model. In order to form a more reliable local model, the RANSAC algorithm is used to filter the points used to fit the model.Although the invalid point does not have the correct depth value, it contains the coordinate informationu,v of the pixel plane. Combining with the local model z^=ax+by+c that has been fitted, the spatial coordinate x,y,z^ of the invalid point can be recalculated, as shown in equation (1). Finally, the invalid point in the depth image is repaired.
(1)x=u−cxz⌢fx,y=v−cyz⌢fy,z⌢=ax+by+c,where cx,cy is the optical center offset of the camera and fx,fy is the focal length of the camera.


### 3.2. Enhancement of the Attitude Calculation

The MARG sensor is a microelectromechanical system (MEMS) composed of a three-axis accelerometer, a three-axis gyroscope, and a three-axis magnetometer. However, due to the impact of sensor manufacturing accuracy and operating environment interference, the raw measurement data of the MARG sensor is generally noisy. From the analysis of its working principle, the gyroscope has good dynamic response characteristics, but it is easy to produce cumulative errors in the attitude calculation, and the accelerometer and magnetometer are just the opposite. Therefore, the QEKF algorithm is used to fuse the three measurement units of the MARG sensor so that the output attitude angle data can be enhanced, as shown in Figure4.Figure 4
Attitude calculation enhancement.
(a)
Block diagram of attitude calculation(b)
The coordinate system of MARG sensorThere are three coordinate systems in this method. Firstly, the acceleration vectoraM=axMayMazMT measured by the accelerometer and the angular velocity vector ωM=ωxMωyMωzMT measured by the gyroscope are both in the body coordinate system. Secondly, the magnetic field intensity vector hM=hxMhyMhzMT measured by the magnetometer is in the coordinate system that can be transformed to the body coordinate system by the rotation matrix Rmb. Lastly, the geomagnetic field intensity vector HG=HxGHyGHzGT and the gravitational acceleration vector g=00−gTG are in the navigation coordinate system that can be transformed from the body coordinate system by the rotation matrix Rbn.Suppose the system state variablesxMk at the k-th time are
(2)xMk=qM0,kqM1,kqM2,kqM3,kbMx,kbMy,kbMz,kT,where q0,kM,q1,kM,q2,kM,q3,kM is the quaternion value of the rotation transformation matrix Rbn at the k-th time and bx,kM,by,kM,bz,kM is the accumulated drift value of the gyroscope on the X, Y, and Z axes at the k-th time.Suppose the system control variablesuMk at the k-th time is
(3)uMk=ωMx,k−bMx,k−1ωMy,k−bMy,k−1ωMz,k−bMz,k−1T.Suppose the system observation variableszMk at the k-th time is
(4)zMk=aMkhMkT.According to the discrete-time model of the quaternion attitude differential equation [23], the following equation can be obtained:
(5)qMk=I4×4+ΩMk−1T2qMk−1,(6)ΩMk−1=0bMx,k−1−ωMx,kbMy,k−1−ωMy,kbMz,k−1−ωMz,kωMx,k−bMx,k−10ωMz,k−bMz,k−1bMy,k−1−ωMy,kωMy,k−bMy,k−1bMz,k−1−ωMz,k0ωMx,k−bMx,k−1ωMz,k−bMz,k−1ωMy,k−bMy,k−1bMx,k−1−ωMx,k0,whereT is the discrete-time interval, Ii×j is the identity matrix with i rows and j columns, qk=q0,kMq1,kMq2,kMq3,kMTM is the quaternion vector of Rbn at the k-th time, and ΩkM is the component of angular velocity in the body coordinate system.Substituting equations (2) and (3) into equation (5), (6) can obtain the state prediction equation:
(7)fMxMk−1,uMk−1=I4×4+ΩMk−1T204×303×4I3×3qMkbMx,k−1bMy,k−1bMz,k−1,where 0i×j is the zero matrix with i rows and j columns.To transform the measured value of the sensor, the following equations can be obtained:(8)aM=Rnb−1gG,(9)hM=Rnb−1Rbm−1HG.Then, the system observation equationhx^kMM was obtained:
(10)a⌢Mx,k=2gq⌢M1,kq⌢M3,k−q⌢M0,kq⌢M2,k,a⌢My,k=2gq⌢M2,kq⌢M3,k+q⌢M0,kq⌢M1,k,a⌢Mz,k=gq⌢M0,k2−q⌢M1,k2−q⌢M2,k2+q⌢M3,k2,h⌢Mx,k=HGxq⌢M0,k2−q⌢M1,k2+q⌢M2,k2−q⌢M3,k2+2HGyq⌢M1,kq⌢M2,k−q⌢M0,kq⌢M3,k−2HGzq⌢M2,kq⌢M3,k+q⌢M0,kq⌢M1,k,h⌢My,k=2HGxq⌢M1,kq⌢M2,k+q⌢M0,kq⌢M3,k+HGyq⌢M0,k2+q⌢M1,k2−q⌢M2,k2−q⌢M3,k2−2HGzq⌢M1,kq⌢M3,k−q⌢M0,kq⌢M2,k,h⌢Mz,k=−2HGxq⌢M2,kq⌢M3,k−q⌢M0,kqM1,k−2HGyq⌢M1,kq⌢M3,k+q⌢M0,kq⌢M2,k+HGzq⌢M0,k2−q⌢M1,k2−q⌢M2,k2+q⌢M3,k2,where aa^x,k,a^y,kM,a^z,kM,h^x,kM,h^y,kM,h^z,kMM is the theoretical measurement value of the sensors and q^0,kM,q^1,kM,q^2,kM,q^3,kM is the quaternion attitude predicted by equation (7).According to the EKF algorithm, calculate the enhanced quaternion attitude at thek-th time:
(11)x⌢Mk=fMxMk−1,uMk−1,P⌢Mk=FMk−1PMk−1FMk−1T+QMk−1,GMk=P⌢MkHMkTHMkP⌢MkHMkT+RMk−1,xMk=x⌢Mk+GMkzMk−hMx⌢Mk,PMk=I−GMkHMkP⌢Mk,where x^kM is the prediction of the system state, P^kM is the prediction of the estimation error covariance matrix PkM,Fk−1M=∂fMxk−1,uk−1/∂xk−1xk−1=xk−1M,HkM=∂hMxk/∂xkxk=x^k−1M, GkM is the Kalman gain, I is the identity matrix, and Qk−1M and RkM are the noise matrix.To further improve the robustness of the system, the process noise matrixQk−1Mand the observation noise matrix RkM in equation (11) adopt an adaptive form.Suppose that the angular velocity measurement valueωM obtained by the gyroscope obeys the Gaussian distribution with variance σg2, and the drift value bM=bxMbyMbzMT obeys the Gaussian distribution with variance σb2, substituting equation (7) to obtain the following:
(12)QMk−1=I7×7σg2+σb2qM1,k−12+qM2,k−12+qM3,k−12T24σg2+σb2qM0,k−12+qM2,k−12+qM3,k−12T24σg2+σb2qM0,k−12+qM1,k−12+qM3,k−12T24σg2+σb2qM0,k−12+qM1,k−12+qM2,k−12T24σb2σb2σb2.Suppose that the acceleration measurement valueaM obtained by the accelerometer obeys the Gaussian distribution with variance σa2, and the magnetic field intensity measurement value hM obtained by the magnetometer obeys the Gaussian distribution with variance σm2. Considering the influence of external acceleration and external magnetic field strength,
(13)RMk=kaaMk−gGk+σa2I3×303×303×3kmhMk−HGk+σm2I3×3,where ka and km are correction coefficients.


## 3.1. Enhancement of the Depth Image

The structured-light RGB-D camera is composed of a color camera, an infrared camera, and an infrared transmitter and obtains the depth value through the speckle feature of the infrared structured light projected on the surface of the object [22]. In a closed and narrow environment, infrared light is easily reflected and absorbed by objects multiple times, so it is easier to produce invalid depth values locally, which in turn introduces additional noise for pose estimation.Since the invalid depth value occurs locally, the invalid point can be repaired by a valid point that is in a similar region to the invalid point. Since the color camera is not affected by the narrow environment, it can be used to determine whether the valid point and the invalid point belong to the same type of patch gathered by the point cloud, that is, the similar region. Therefore, the enhanced method proposed in this paper includes two steps, namely, similar region extraction and local point cloud reconstruction, as shown in Figure3. The former extracts similar regions in the depth image through the texture features of the registered color image. The latter uses valid points in the similar region to construct a local model of the region and then recalculates invalid points in the region with the model to complete the enhancement repair process.Figure 3
Depth image enhancement in a narrow environment.In the similar region extraction step, the Mean Shift algorithm is used to extract the color texture information of the color image, and the Flood Fill algorithm is used to extract the corresponding color-connected domains to form a similar region mask. Finally, the mask is applied to the depth image to obtain the region segmentation of the similar point cloud. Since the Mean Shift algorithm needs to perform multiple iterative calculations for each pixel, it will affect the real-time performance of the system in practical applications. In order to speed up the calculation process, this paper downsamples the input color image to a lower resolution, then forms the similar region mask on the low-resolution image, and finally performs bilinear interpolation on the mask to restore the original scale.In the local point cloud reconstruction step, by fitting the patch model in the similar region and recalculating the depth value of the invalid point on the patch, the invalid point is repaired. The least-squares method is the basic method to fit the patch model, but in practical applications, the region segmented by the color image has deviation, and not all valid points can be used to fit the patch model. In order to form a more reliable local model, the RANSAC algorithm is used to filter the points used to fit the model.Although the invalid point does not have the correct depth value, it contains the coordinate informationu,v of the pixel plane. Combining with the local model z^=ax+by+c that has been fitted, the spatial coordinate x,y,z^ of the invalid point can be recalculated, as shown in equation (1). Finally, the invalid point in the depth image is repaired.
(1)x=u−cxz⌢fx,y=v−cyz⌢fy,z⌢=ax+by+c,where cx,cy is the optical center offset of the camera and fx,fy is the focal length of the camera.


## 3.2. Enhancement of the Attitude Calculation

The MARG sensor is a microelectromechanical system (MEMS) composed of a three-axis accelerometer, a three-axis gyroscope, and a three-axis magnetometer. However, due to the impact of sensor manufacturing accuracy and operating environment interference, the raw measurement data of the MARG sensor is generally noisy. From the analysis of its working principle, the gyroscope has good dynamic response characteristics, but it is easy to produce cumulative errors in the attitude calculation, and the accelerometer and magnetometer are just the opposite. Therefore, the QEKF algorithm is used to fuse the three measurement units of the MARG sensor so that the output attitude angle data can be enhanced, as shown in Figure4.Figure 4
Attitude calculation enhancement.
(a)
Block diagram of attitude calculation(b)
The coordinate system of MARG sensorThere are three coordinate systems in this method. Firstly, the acceleration vectoraM=axMayMazMT measured by the accelerometer and the angular velocity vector ωM=ωxMωyMωzMT measured by the gyroscope are both in the body coordinate system. Secondly, the magnetic field intensity vector hM=hxMhyMhzMT measured by the magnetometer is in the coordinate system that can be transformed to the body coordinate system by the rotation matrix Rmb. Lastly, the geomagnetic field intensity vector HG=HxGHyGHzGT and the gravitational acceleration vector g=00−gTG are in the navigation coordinate system that can be transformed from the body coordinate system by the rotation matrix Rbn.Suppose the system state variablesxMk at the k-th time are
(2)xMk=qM0,kqM1,kqM2,kqM3,kbMx,kbMy,kbMz,kT,where q0,kM,q1,kM,q2,kM,q3,kM is the quaternion value of the rotation transformation matrix Rbn at the k-th time and bx,kM,by,kM,bz,kM is the accumulated drift value of the gyroscope on the X, Y, and Z axes at the k-th time.Suppose the system control variablesuMk at the k-th time is
(3)uMk=ωMx,k−bMx,k−1ωMy,k−bMy,k−1ωMz,k−bMz,k−1T.Suppose the system observation variableszMk at the k-th time is
(4)zMk=aMkhMkT.According to the discrete-time model of the quaternion attitude differential equation [23], the following equation can be obtained:
(5)qMk=I4×4+ΩMk−1T2qMk−1,(6)ΩMk−1=0bMx,k−1−ωMx,kbMy,k−1−ωMy,kbMz,k−1−ωMz,kωMx,k−bMx,k−10ωMz,k−bMz,k−1bMy,k−1−ωMy,kωMy,k−bMy,k−1bMz,k−1−ωMz,k0ωMx,k−bMx,k−1ωMz,k−bMz,k−1ωMy,k−bMy,k−1bMx,k−1−ωMx,k0,whereT is the discrete-time interval, Ii×j is the identity matrix with i rows and j columns, qk=q0,kMq1,kMq2,kMq3,kMTM is the quaternion vector of Rbn at the k-th time, and ΩkM is the component of angular velocity in the body coordinate system.Substituting equations (2) and (3) into equation (5), (6) can obtain the state prediction equation:
(7)fMxMk−1,uMk−1=I4×4+ΩMk−1T204×303×4I3×3qMkbMx,k−1bMy,k−1bMz,k−1,where 0i×j is the zero matrix with i rows and j columns.To transform the measured value of the sensor, the following equations can be obtained:(8)aM=Rnb−1gG,(9)hM=Rnb−1Rbm−1HG.Then, the system observation equationhx^kMM was obtained:
(10)a⌢Mx,k=2gq⌢M1,kq⌢M3,k−q⌢M0,kq⌢M2,k,a⌢My,k=2gq⌢M2,kq⌢M3,k+q⌢M0,kq⌢M1,k,a⌢Mz,k=gq⌢M0,k2−q⌢M1,k2−q⌢M2,k2+q⌢M3,k2,h⌢Mx,k=HGxq⌢M0,k2−q⌢M1,k2+q⌢M2,k2−q⌢M3,k2+2HGyq⌢M1,kq⌢M2,k−q⌢M0,kq⌢M3,k−2HGzq⌢M2,kq⌢M3,k+q⌢M0,kq⌢M1,k,h⌢My,k=2HGxq⌢M1,kq⌢M2,k+q⌢M0,kq⌢M3,k+HGyq⌢M0,k2+q⌢M1,k2−q⌢M2,k2−q⌢M3,k2−2HGzq⌢M1,kq⌢M3,k−q⌢M0,kq⌢M2,k,h⌢Mz,k=−2HGxq⌢M2,kq⌢M3,k−q⌢M0,kqM1,k−2HGyq⌢M1,kq⌢M3,k+q⌢M0,kq⌢M2,k+HGzq⌢M0,k2−q⌢M1,k2−q⌢M2,k2+q⌢M3,k2,where aa^x,k,a^y,kM,a^z,kM,h^x,kM,h^y,kM,h^z,kMM is the theoretical measurement value of the sensors and q^0,kM,q^1,kM,q^2,kM,q^3,kM is the quaternion attitude predicted by equation (7).According to the EKF algorithm, calculate the enhanced quaternion attitude at thek-th time:
(11)x⌢Mk=fMxMk−1,uMk−1,P⌢Mk=FMk−1PMk−1FMk−1T+QMk−1,GMk=P⌢MkHMkTHMkP⌢MkHMkT+RMk−1,xMk=x⌢Mk+GMkzMk−hMx⌢Mk,PMk=I−GMkHMkP⌢Mk,where x^kM is the prediction of the system state, P^kM is the prediction of the estimation error covariance matrix PkM,Fk−1M=∂fMxk−1,uk−1/∂xk−1xk−1=xk−1M,HkM=∂hMxk/∂xkxk=x^k−1M, GkM is the Kalman gain, I is the identity matrix, and Qk−1M and RkM are the noise matrix.To further improve the robustness of the system, the process noise matrixQk−1Mand the observation noise matrix RkM in equation (11) adopt an adaptive form.Suppose that the angular velocity measurement valueωM obtained by the gyroscope obeys the Gaussian distribution with variance σg2, and the drift value bM=bxMbyMbzMT obeys the Gaussian distribution with variance σb2, substituting equation (7) to obtain the following:
(12)QMk−1=I7×7σg2+σb2qM1,k−12+qM2,k−12+qM3,k−12T24σg2+σb2qM0,k−12+qM2,k−12+qM3,k−12T24σg2+σb2qM0,k−12+qM1,k−12+qM3,k−12T24σg2+σb2qM0,k−12+qM1,k−12+qM2,k−12T24σb2σb2σb2.Suppose that the acceleration measurement valueaM obtained by the accelerometer obeys the Gaussian distribution with variance σa2, and the magnetic field intensity measurement value hM obtained by the magnetometer obeys the Gaussian distribution with variance σm2. Considering the influence of external acceleration and external magnetic field strength,
(13)RMk=kaaMk−gGk+σa2I3×303×303×3kmhMk−HGk+σm2I3×3,where ka and km are correction coefficients.


## 4. Sensor Fusion Based on Modified RBF-AUKF Method

### 4.1. End-Effector Pose Prediction Model

The RGB-D camera and MARG sensor are both arranged at the end of the manipulator, so the sensor data after preenhancement processing is directly in the end-effector coordinate system. The joint encoders are arranged at each joint of the manipulator, so it is necessary to transform the angle of each joint to the pose in the end-effector coordinate system through forward kinematics calculation.Suppose the system state variablesxk at the k-th time is
(14)xk=pkqkT,where pk=px,kpy,kpz,kT is the position of the end-effector at the k-th time and qk=q0,kq1,kq2,kq3,kT is the quaternion attitude of the end-effector at the k-th time.Suppose the system control variablesuk at the k-th time is
(15)uk=ωkṗkp¨kT,where ωk=ωx,kωy,kωz,kT is the X, Y, and Z axis angular velocities of the end-effector, ṗk=ṗx,kṗy,kṗz,kT is the position velocity of the end-effector at the k-th time, and p¨k=p¨x,kp¨y,kp¨z,kT is the position acceleration of the end-effector at the k-th time.According to the differential model of motion [24], the system state prediction equation fxk−1,uk−1 can be obtained:
(16)pk=pk−1+Rnb,k−1ṗx,k−1T+p¨x,k−1T22ṗy,k−1T+p¨y,k−1T22ṗz,k−1T+p¨z,k−1T22,qk=qk−1+T20−ωx,k−1−ωy,k−1−ωz,k−1ωx,k−10ωz,k−1−ωy,k−1ωy,k−1−ωz,k−10ωx,k−1ωz,k−1ωy,k−1−ωx,k−10qk−1,where T is the discrete-time interval and Rb,k−1n is the rotation matrix from the body coordinate system to the navigation coordinate system as shown in equation (17).
(17)Rnb,k−1=q0,k−12+q1,k−12−q2,k−12−q3,k−122q1,k−1q2,k−1−q0,k−1q3,k−12q1,k−1q3,k−1+q0,k−1q2,k−12q1,k−1q2,k−1+q0,k−1q3,k−1q0,k−12−q1,k−12+q2,k−12−q3,k−122q2,k−1q3,k−1−q0,k−1q1,k−12q1,k−1q3,k−1−q0,k−1q2,k−12q2,k−1q3,k−1+q0,k−1q1,k−1q0,k−12−q1,k−12−q2,k−12+q3,k−12.The state prediction equation mentioned in equation (16) has obvious nonlinear characteristics. Both EKF and UKF algorithms can solve the problem of nonlinear state estimation. However, when the calculation performance can meet the requirements, the estimation accuracy of the UKF algorithm is higher [25]. By sampling, the UKF algorithm approximates the probability density distribution of the nonlinear prediction equation, and the constructed sampling point set is as follows:
(18)χki=xk,i=0,xk+N+λPk−1i,i=1,⋯,N,xk−N+λPk−1i,i=N+1,⋯,2N,where χki is the i-th column of the sample point set matrix at the k-th time, which is a 7-dimensional column vector, N+λPk−1i is the i-th column of the matrix N+λPk−1, which is an N-dimensional column vector, λ is the scale parameter to adjust the approximation accuracy, and Pk−1 is the estimation error covariance matrix.


### 4.2. Multirate Update Model

Due to the limitation of sensor processing capacity, the measurement update frequency of different sensors is not the same, or even quite different, which means that even if the sensors are triggered to collect at the same time, the measurement of each sensor cannot be received at the same time.Although it is difficult for sensors to publish data synchronously all the time, it is possible to approximate the software synchronization of multiple sensors through timestamp alignment. The basic strategy is to set a sliding time window; when all sensor data falls within this window, it can be considered that the sensor data is synchronous. As is shown in Figure5, d is the size of the sliding time window, which determines the accuracy of data synchronization. Its maximum size is the smallest update time interval among all sensors. The smaller the size, the higher the accuracy of the time, but the lower the frequency of synchronous data.Figure 5
An update model that maintains an average update frequency.The UKF algorithm corrects the prediction model through sensor observations. Therefore, when the algorithm is used for multisensor fusion, the fusion accuracy is greatly affected by the time synchronization [26]. Because the update frequency of synchronous data is generally much lower than the average update frequency of data, the real-time performance of synchronization fusion estimation is poor. To solve this shortcoming, a modified synchronization fusion algorithm is proposed including two update models, namely, the absolute state update model and the relative state update model.Define a matrixHm×n with m rows and n columns as a state update matrix, and its i-th row and j-th column satisfy the following:
(19)Hm×ni,j=1,measurementiisrelatedtostatej,0,otherwise.The absolute state update model updates the system state variables including the position and attitude. The number of rows and columns in the state update matrixHm×n satisfies the following:
(20)m=∑k=1NAAsk,n=7,where NA is the total number of sensors synchronized with the system state variables and NA must be equal to Nmax which is the maximum number of sensors; otherwise Hm×n=0. skA is the total number of the system state variables that the k-th sensor can observe.The relative state update model updates the system control variables including velocity and acceleration. The number of rows and columns in the state update matrixHm×nsatisfies the following:
(21)m=∑k=1NRRsk,n=9,where NR is the total number of sensors synchronized with the system control variables and NR must satisfy 2≤NR≤Nmax; otherwise Hm×n=0. skR is the total number of system control variables that the k-th sensor can observe.When the multisensor data is not all synchronized, the relative state update is used to increase the update frequency locally, and when the data is all synchronized, the system state is globally corrected to improve the estimation accuracy. The final observation equation is as follows:(22)hxk∣k−1,uk=HA00HRxk∣k−1uk,where HA is the absolute state update matrix determined by equation (20) and HR is the relative state update matrix determined by equation (21).According to the UKF algorithm, calculate the fusion pose estimationxk of the end-effector:
(23)χk∣k−1=fχk−1,uk−1,xk∣k−1=∑i=02nWiχk∣k−1i,Px,k∣k−1=∑i=02nWiχk∣k−1i−xk∣k−1χk∣k−1i−xk∣k−1T+Qk,γk∣k−1=hχk∣k−1,uk,zk∣k−1=∑i=02nWiγk∣k−1i,Pz,k∣k−1=∑i=02nWiγk∣k−1i−zk∣k−1γk∣k−1i−zk∣k−1T+Rk,Pxz,k∣k−1=∑i=02nWiχk∣k−1i−xk∣k−1γk∣k−1i−zk∣k−1T,Gk=Pxz,k∣k−1Pz,k∣k−1−1,xk=xk∣k−1+Gkzk−zk∣k−1,Pk=Px,k∣k−1−GkPz,k∣k−1GkT,where χk∣k−1 is the state prediction of sampling points determined by equation (16), xk∣k−1 is the weighted mean of the state prediction, Wi is the weight value determined by equation (24), Px,k∣k−1 is the weighted variance of the state prediction, γk∣k−1 is the state observation of sampling points determined by equation (22), zk∣k−1 is the weighted mean of the state observation, Pz,k∣k−1 is the weighted variance of the state observation, Pxz,k∣k−1 is the weighted variance of the state prediction and observation, zk is the sensor measurement value, Gk is the Kalman gain, Qk and Rk are the noise matrix.
(24)Wi=NN+λ,i=0N,N2N+λ,i=1,⋯,2N.


### 4.3. Adaptive Adjustment of the Noise Matrix

It can be seen from equation (23) that the process noise covariance matrix Qk and the observation noise covariance matrix Rk can affect the confidence level of the data during state prediction and update. In general applications, the variance of system state noise and sensor measurement noise is a fixed value that can be counted in advance. However, in a complex environment, the variance of noise will fluctuate to a certain extent. If the fixed value is still set in advance, it will affect the accuracy of the final fusion state estimation. Therefore, the system identification ability of the RBF network needs to be used to adaptively adjust the noise variance.The RBF network includes a three-layer structure of input layer, hidden layer, and output layer, as shown in Figure6. The prediction error xk−xk∣k−1, estimation error covariance matrix Pk, Kalman gain Gk, and observation error zk−zk∣k−1 shown in equation (23) can reflect the deviation and fluctuation during the prediction and update process in the UKF and can be used as features of the input layer.Figure 6
RBF network for adjusting noise matrix.Considering the similarity between the features, select the features according to the position state, attitude state, and the observations of each sensor:(25)xk−xk∣k−1=ep,keq,kT,Pk=Ppp,kPpq,kPqp,kPqq,k,Gk=Gs1,kGs2,kGs3,k,zk−zk∣k−1=es1,kes2,kes3,kT,where ep,k,eq,k are the position and attitude prediction error; Ppp,k,Ppq,k,Pqp,k,Pqq,k are the estimation error covariance matrix of position, position-attitude, attitude-position, and attitude; Gs1,k,Gs2,k,Gs3,k are the Kalman gain of the three sensors; and es1,k,es2,k,es3,k are the observation error of the three sensors.Extract the main information of the feature matrix in equation (25) to obtain the input layer feature vector xinput, which contains 10 feature information:
(26)xinput=ep,keq,ktrPpp,ktrPqq,kGs1,kGs2,kGs3,kes1,kes2,kes3,kT.From input layer to hidden layer, radial basis function is used to activate neurons:(27)ϕixinput=exinput−μi/σi2,where μi,σi are the center and width of the i-th hidden layer neuron.The linear function is used from hidden layer to output layer:(28)youtput,j=∑iωijϕixinput,where ωij is the weight of the i-th hidden layer neuron to the j-th output layer neuron.The output of the output layer is as follows:(29)youtput=σpσqσs1σs2σs3T,(30)Qk+1=IσpσqT,Rk+1=Iσs1σs2σs3T,where σp and σq are the standard deviations of the process noise of the position and attitude and σs1, σs2, and σs3 are the standard deviations of the observation noise of the three sensors.The training process of the network is divided into two stages: offline and online. The offline training stage uses the data with known noise variance to obtain the initial neuron center and connection weights. At this stage, the output error of the network is the difference between the true noise variance and the network predicted noise variance. Since the noise variance in the offline stage is not the noise variance in the actual experiment, it is easy to produce overfitting in the offline training stage, which needs to be corrected through the online stage. The online training stage uses the updated prediction errorxk+1−xk+1∣k at the next time as the output error of the network. The training goal is to minimize the prediction error:
(31)Loss=argminμi,ωijxk+1−xk+1∣k.


## 4.1. End-Effector Pose Prediction Model

The RGB-D camera and MARG sensor are both arranged at the end of the manipulator, so the sensor data after preenhancement processing is directly in the end-effector coordinate system. The joint encoders are arranged at each joint of the manipulator, so it is necessary to transform the angle of each joint to the pose in the end-effector coordinate system through forward kinematics calculation.Suppose the system state variablesxk at the k-th time is
(14)xk=pkqkT,where pk=px,kpy,kpz,kT is the position of the end-effector at the k-th time and qk=q0,kq1,kq2,kq3,kT is the quaternion attitude of the end-effector at the k-th time.Suppose the system control variablesuk at the k-th time is
(15)uk=ωkṗkp¨kT,where ωk=ωx,kωy,kωz,kT is the X, Y, and Z axis angular velocities of the end-effector, ṗk=ṗx,kṗy,kṗz,kT is the position velocity of the end-effector at the k-th time, and p¨k=p¨x,kp¨y,kp¨z,kT is the position acceleration of the end-effector at the k-th time.According to the differential model of motion [24], the system state prediction equation fxk−1,uk−1 can be obtained:
(16)pk=pk−1+Rnb,k−1ṗx,k−1T+p¨x,k−1T22ṗy,k−1T+p¨y,k−1T22ṗz,k−1T+p¨z,k−1T22,qk=qk−1+T20−ωx,k−1−ωy,k−1−ωz,k−1ωx,k−10ωz,k−1−ωy,k−1ωy,k−1−ωz,k−10ωx,k−1ωz,k−1ωy,k−1−ωx,k−10qk−1,where T is the discrete-time interval and Rb,k−1n is the rotation matrix from the body coordinate system to the navigation coordinate system as shown in equation (17).
(17)Rnb,k−1=q0,k−12+q1,k−12−q2,k−12−q3,k−122q1,k−1q2,k−1−q0,k−1q3,k−12q1,k−1q3,k−1+q0,k−1q2,k−12q1,k−1q2,k−1+q0,k−1q3,k−1q0,k−12−q1,k−12+q2,k−12−q3,k−122q2,k−1q3,k−1−q0,k−1q1,k−12q1,k−1q3,k−1−q0,k−1q2,k−12q2,k−1q3,k−1+q0,k−1q1,k−1q0,k−12−q1,k−12−q2,k−12+q3,k−12.The state prediction equation mentioned in equation (16) has obvious nonlinear characteristics. Both EKF and UKF algorithms can solve the problem of nonlinear state estimation. However, when the calculation performance can meet the requirements, the estimation accuracy of the UKF algorithm is higher [25]. By sampling, the UKF algorithm approximates the probability density distribution of the nonlinear prediction equation, and the constructed sampling point set is as follows:
(18)χki=xk,i=0,xk+N+λPk−1i,i=1,⋯,N,xk−N+λPk−1i,i=N+1,⋯,2N,where χki is the i-th column of the sample point set matrix at the k-th time, which is a 7-dimensional column vector, N+λPk−1i is the i-th column of the matrix N+λPk−1, which is an N-dimensional column vector, λ is the scale parameter to adjust the approximation accuracy, and Pk−1 is the estimation error covariance matrix.


## 4.2. Multirate Update Model

Due to the limitation of sensor processing capacity, the measurement update frequency of different sensors is not the same, or even quite different, which means that even if the sensors are triggered to collect at the same time, the measurement of each sensor cannot be received at the same time.Although it is difficult for sensors to publish data synchronously all the time, it is possible to approximate the software synchronization of multiple sensors through timestamp alignment. The basic strategy is to set a sliding time window; when all sensor data falls within this window, it can be considered that the sensor data is synchronous. As is shown in Figure5, d is the size of the sliding time window, which determines the accuracy of data synchronization. Its maximum size is the smallest update time interval among all sensors. The smaller the size, the higher the accuracy of the time, but the lower the frequency of synchronous data.Figure 5
An update model that maintains an average update frequency.The UKF algorithm corrects the prediction model through sensor observations. Therefore, when the algorithm is used for multisensor fusion, the fusion accuracy is greatly affected by the time synchronization [26]. Because the update frequency of synchronous data is generally much lower than the average update frequency of data, the real-time performance of synchronization fusion estimation is poor. To solve this shortcoming, a modified synchronization fusion algorithm is proposed including two update models, namely, the absolute state update model and the relative state update model.Define a matrixHm×n with m rows and n columns as a state update matrix, and its i-th row and j-th column satisfy the following:
(19)Hm×ni,j=1,measurementiisrelatedtostatej,0,otherwise.The absolute state update model updates the system state variables including the position and attitude. The number of rows and columns in the state update matrixHm×n satisfies the following:
(20)m=∑k=1NAAsk,n=7,where NA is the total number of sensors synchronized with the system state variables and NA must be equal to Nmax which is the maximum number of sensors; otherwise Hm×n=0. skA is the total number of the system state variables that the k-th sensor can observe.The relative state update model updates the system control variables including velocity and acceleration. The number of rows and columns in the state update matrixHm×nsatisfies the following:
(21)m=∑k=1NRRsk,n=9,where NR is the total number of sensors synchronized with the system control variables and NR must satisfy 2≤NR≤Nmax; otherwise Hm×n=0. skR is the total number of system control variables that the k-th sensor can observe.When the multisensor data is not all synchronized, the relative state update is used to increase the update frequency locally, and when the data is all synchronized, the system state is globally corrected to improve the estimation accuracy. The final observation equation is as follows:(22)hxk∣k−1,uk=HA00HRxk∣k−1uk,where HA is the absolute state update matrix determined by equation (20) and HR is the relative state update matrix determined by equation (21).According to the UKF algorithm, calculate the fusion pose estimationxk of the end-effector:
(23)χk∣k−1=fχk−1,uk−1,xk∣k−1=∑i=02nWiχk∣k−1i,Px,k∣k−1=∑i=02nWiχk∣k−1i−xk∣k−1χk∣k−1i−xk∣k−1T+Qk,γk∣k−1=hχk∣k−1,uk,zk∣k−1=∑i=02nWiγk∣k−1i,Pz,k∣k−1=∑i=02nWiγk∣k−1i−zk∣k−1γk∣k−1i−zk∣k−1T+Rk,Pxz,k∣k−1=∑i=02nWiχk∣k−1i−xk∣k−1γk∣k−1i−zk∣k−1T,Gk=Pxz,k∣k−1Pz,k∣k−1−1,xk=xk∣k−1+Gkzk−zk∣k−1,Pk=Px,k∣k−1−GkPz,k∣k−1GkT,where χk∣k−1 is the state prediction of sampling points determined by equation (16), xk∣k−1 is the weighted mean of the state prediction, Wi is the weight value determined by equation (24), Px,k∣k−1 is the weighted variance of the state prediction, γk∣k−1 is the state observation of sampling points determined by equation (22), zk∣k−1 is the weighted mean of the state observation, Pz,k∣k−1 is the weighted variance of the state observation, Pxz,k∣k−1 is the weighted variance of the state prediction and observation, zk is the sensor measurement value, Gk is the Kalman gain, Qk and Rk are the noise matrix.
(24)Wi=NN+λ,i=0N,N2N+λ,i=1,⋯,2N.


## 4.3. Adaptive Adjustment of the Noise Matrix

It can be seen from equation (23) that the process noise covariance matrix Qk and the observation noise covariance matrix Rk can affect the confidence level of the data during state prediction and update. In general applications, the variance of system state noise and sensor measurement noise is a fixed value that can be counted in advance. However, in a complex environment, the variance of noise will fluctuate to a certain extent. If the fixed value is still set in advance, it will affect the accuracy of the final fusion state estimation. Therefore, the system identification ability of the RBF network needs to be used to adaptively adjust the noise variance.The RBF network includes a three-layer structure of input layer, hidden layer, and output layer, as shown in Figure6. The prediction error xk−xk∣k−1, estimation error covariance matrix Pk, Kalman gain Gk, and observation error zk−zk∣k−1 shown in equation (23) can reflect the deviation and fluctuation during the prediction and update process in the UKF and can be used as features of the input layer.Figure 6
RBF network for adjusting noise matrix.Considering the similarity between the features, select the features according to the position state, attitude state, and the observations of each sensor:(25)xk−xk∣k−1=ep,keq,kT,Pk=Ppp,kPpq,kPqp,kPqq,k,Gk=Gs1,kGs2,kGs3,k,zk−zk∣k−1=es1,kes2,kes3,kT,where ep,k,eq,k are the position and attitude prediction error; Ppp,k,Ppq,k,Pqp,k,Pqq,k are the estimation error covariance matrix of position, position-attitude, attitude-position, and attitude; Gs1,k,Gs2,k,Gs3,k are the Kalman gain of the three sensors; and es1,k,es2,k,es3,k are the observation error of the three sensors.Extract the main information of the feature matrix in equation (25) to obtain the input layer feature vector xinput, which contains 10 feature information:
(26)xinput=ep,keq,ktrPpp,ktrPqq,kGs1,kGs2,kGs3,kes1,kes2,kes3,kT.From input layer to hidden layer, radial basis function is used to activate neurons:(27)ϕixinput=exinput−μi/σi2,where μi,σi are the center and width of the i-th hidden layer neuron.The linear function is used from hidden layer to output layer:(28)youtput,j=∑iωijϕixinput,where ωij is the weight of the i-th hidden layer neuron to the j-th output layer neuron.The output of the output layer is as follows:(29)youtput=σpσqσs1σs2σs3T,(30)Qk+1=IσpσqT,Rk+1=Iσs1σs2σs3T,where σp and σq are the standard deviations of the process noise of the position and attitude and σs1, σs2, and σs3 are the standard deviations of the observation noise of the three sensors.The training process of the network is divided into two stages: offline and online. The offline training stage uses the data with known noise variance to obtain the initial neuron center and connection weights. At this stage, the output error of the network is the difference between the true noise variance and the network predicted noise variance. Since the noise variance in the offline stage is not the noise variance in the actual experiment, it is easy to produce overfitting in the offline training stage, which needs to be corrected through the online stage. The online training stage uses the updated prediction errorxk+1−xk+1∣k at the next time as the output error of the network. The training goal is to minimize the prediction error:
(31)Loss=argminμi,ωijxk+1−xk+1∣k.


## 5. Experiments and Results

### 5.1. Experimental Device

The mechanical structure of the experimental device consisted of a 1-DOF motion platform and an 8-DOF redundant manipulator, as shown in Figure7. The joint encoders were arranged at each rotary joint, and the RGB-D camera and the MARG sensor were arranged at the end of the manipulator. The RGB-D camera was connected to the industrial personal computer (IPC) through USB, and the MARG sensor and the joint encoders were connected to the sensor collector through the 485 bus and then connected to the IPC through USB.Figure 7
Experimental device.


### 5.2. Experiment of the Enhancement Part

In order to test the effect of the sensor complementary enhancement method proposed in this article, the experimental device was placed in a narrow artificially built environment, considering that the RGB-D camera focuses more on static repeatability and the MARG sensor focuses more on dynamic stability. Therefore, when testing the enhancement effect of the RGB-D camera, a periodically repeated expected motion trajectory was used, while the MARG sensor used a random expected trajectory. The experimental results are shown in Figures8 and 9.Figure 8
Results of depth image enhancement.
(a)(b)(c)Figure 9
Results of attitude calculation enhancement.Comparing Figures8(a) and 8(b), it can be seen that the point cloud after enhancement is closer to the real environment and has a less invalid blind area caused by the narrow environment. Figure 8(c) further demonstrates the impact of depth image enhancement on the accuracy of position estimation. It can be seen from the figure that the position trajectory obtained by using the original image information fluctuates more drastically and deviates from the reference trajectory to a greater degree. But the position trajectory obtained by using the enhanced image information is closer to the reference trajectory, and the overall fluctuation is smaller, which is more conducive to the fusion pose estimation with other sensors.Comparing the desired and enhanced quaternion attitude in Figure9, the overall calculation error of the QEKF algorithm proposed in this paper is within the expected range, which means acceleration information can compensate for the integral drift of angular velocity and geomagnetic field information can calibrate angle measurement errors. Through the enhancement of attitude calculation, it can provide more reliable attitude data for multisensor fusion.


### 5.3. Simulation of the Fusion Part

Because the statistical tests are very sensitive when the sample size is small, both simulation tests and actual experiments were carried out to complete the verification of the multisensor fusion pose estimation method. Simulation tests were mainly used to test the stability of the method under extreme conditions, and actual experiments were mainly used to test the fusion effect of the method in real environments. It was easy to adjust the sensor parameters in simulation environments, so a large number of samples can be generated to test the performance of the method. Although the experimental samples that were collected in the real environments were limited, these samples were more targeted and reliable after the screening of the experimental methods by the simulation tests. At the same time, to present the experimental results as objectively as possible, the root mean square (RMS) errors of multiple experimental results were used to evaluate the performance differences of different methods.Since real sensors can only be simulated numerically under simulation conditions, the actual state and sensors were simplified during the simulation. Different from the real system, an observation state and two sensors were set up in the simulation environment. The observation noise of sensor 1 was smaller than that of sensor 2, and sensor 1 only produced a stable zero offset, while the zero offset of sensor 2 increased with time. Figure10 and Table 1 show the effect of the synchronization update method on the fusion result.Figure 10
The impact of synchronization.
(a)
Single sensor 1(b)
Single sensor 2(c)
Sensor 1 and 2 update frequency:f1>f2(d)
Sensor 1 and 2 update frequency:f1<f2(e)
Sensor 1 and 2 synchronizedTable 1
RMS errors of different update methods.
Update methodSensor1SingleSensor2SingleSensors 1 and 2f1>f2Sensors 1 and 2f1>f2Sensors 1 and 2SynchronizedRMS error1.021.160.8240.9620.413Since the noise of sensor 2 is greater than that of sensor 1, the RMS error of sensor 2 is greater than that of sensor 1 when the state of a single sensor is estimated. In the simulation, the data update frequencyf1 of sensor 1 was different from the data update frequency f2 of sensor 2. If the data synchronization of the sensor was not considered, the state update was performed when the data of any sensor inputted. From the experimental results shown in Figures 10(c) and 10(d), the RMS error when f1>f2 is less than that when f1<f2, indicating that the fusion result is affected by the sensor update frequency, and it is biased towards high-sampling-rate data sources. If the high-sampling-rate data is noisy, the fusion accuracy cannot be significantly improved. The above defects are improved in the modified synchronization update method. It can be seen from the experimental results shown in Figure 10(e) that the RMS error of the fusion data is significantly reduced, and there will be no sawtooth fluctuations in the asynchronous fusion. Therefore, the data synchronization of sensors is an important factor in multirate sensor fusion.In the simulation environment, the process and observation noise conformed to the Gaussian distribution, but the variance of the distribution changed at every moment and satisfied a nonlinear model. The RBF network was designed to predict the variances of the process and observation noise and was trained when the filter was updated. In the simulation, the variances predicted by the network in the test set were compared with the preset variances, and the coefficient of determination (marked asR-square) was calculated. The closer the coefficient of determination is to 1, the higher the accuracy of the network. It can be seen from the results shown in Figure 11 that the network can predict the variances of the noise accurately and has no low coefficients of determination (0.71 for the process noise and 0.84 for the observation noise).Figure 11
Prediction of the noise matrix.Although the coefficient of determination cannot be obtained when the true variance is unknown, the accuracy of the network can be indirectly reflected through the state estimation deviation. In other words, if the state estimation deviation can remain stable even in extreme cases, it means that the network’s prediction of the parameters is accurate. In the simulation, the process noise and the observation noise of the sensor were designed to increase with time. Figure12 and Table 2 show the fusion effects of different fusion algorithms in this case.Figure 12
The fusion effects of different fusion algorithms.
(a)(b)Table 2
RMS errors of different fusion methods.
Fusion methodEKFUKFRBF-AEKFRBF-AUKFRMS error39.4633.4214.987.863It can be seen from the experimental results that at the beginning, the system noise is small and the estimation error of each fusion method is not much different, but as time goes by, the system noise increases. The original EKF and UKF algorithms cannot adaptively adjust the filter parameters, which are easily affected by noise and produce state estimation deviations. Through the RBF network proposed in this paper, the noise covariance matrix can be predicted, so that the EKF and UKF algorithms have adaptive adjustment capabilities. The modified RBF-AEKF and RBF-AUKF methods have smaller RMS errors than the original EKF and UKF methods, so they can be more robust in complex environments.


### 5.4. Experiment of the Fusion Part

Since the environment established in the simulation can only approximately reflect the actual situation, it is still necessary to conduct experiments in the real environment. In actual experiments, the end-effector of the manipulator periodically moved around theZ-axis of the base coordinate system, and the pose of the end-effector was set as the desired value. The position estimations px,py,pz along the X, Y, and Z axes were obtained by fusion of the visual odometer based on the RGB-D camera and the kinematics calculation based on joint encoders. And the quaternion attitude estimations q0,q1,q2,q3 were obtained by fusion of the MARG sensor, visual odometer, and encoders. The experiment compares the original UKF method and the modified adaptive UKF method proposed in this paper, and the results are shown in Figures 13 and 14 and Table 3.Figure 13
Results of position estimation.Figure 14
Results of attitude estimation.Table 3
RMS errors of the original method and modified method.
RMS errorEstimated variablepx(cm)py(cm)pz(cm)q0q1q2q3MARG———0.3190.1010.2600.273VO2.767.022.190.1090.1940.2200.232Encoder0.3031.080.3600.09330.06240.05840.133Fusion (modified)0.1750.3640.1080.008050.005640.006450.0482Fusion (original)0.6631.590.4750.008050.005640.006450.0482It can be seen from Table3 that the RMS error of the pose estimation after multisensor fusion is smaller than that of a single sensor. Compared with that of the original method, the RMS error of the modified method proposed in this paper is significantly reduced in position estimation, but there is little difference in attitude estimation. The main reason is the modified method proposed in this paper is better than the traditional method due to its better adaptive effect when the noise fluctuation is large. In this experiment, the overall fluctuation of the sensor in attitude measurement is smaller than that of position measurement, so the modified method is close to the original method in attitude estimation.It can be seen from Figure13 that the visual odometer is prone to local sudden changes due to environmental interference and the encoder is prone to measurement offset. Compared with the original method, the modified method handles the measurement data of local sudden change better, so the overall RMS error is smaller. And because of the combination of absolute and relative updates, the overall update rate is maintained. In this experiment, the original method was updated only 258 times, while the modified method was updated 479 times, so the real-time response rate of state estimation was improved. It can be seen from Figure 14 that the rotation angle of the desired trajectory on the X and Y axes is small, and the rotation angle on the Z-axis is large. When the rotation angle is small, it can be seen that the sensor has obvious measurement drift, when the rotation angle is large, it is easy to produce local sudden changes. The modified method in this paper can compensate the measurement drift of each sensor, reduce the influence of local sudden changes on the state estimation, and thus can better follow the desired trajectory.


## 5.1. Experimental Device

The mechanical structure of the experimental device consisted of a 1-DOF motion platform and an 8-DOF redundant manipulator, as shown in Figure7. The joint encoders were arranged at each rotary joint, and the RGB-D camera and the MARG sensor were arranged at the end of the manipulator. The RGB-D camera was connected to the industrial personal computer (IPC) through USB, and the MARG sensor and the joint encoders were connected to the sensor collector through the 485 bus and then connected to the IPC through USB.Figure 7
Experimental device.


## 5.2. Experiment of the Enhancement Part

In order to test the effect of the sensor complementary enhancement method proposed in this article, the experimental device was placed in a narrow artificially built environment, considering that the RGB-D camera focuses more on static repeatability and the MARG sensor focuses more on dynamic stability. Therefore, when testing the enhancement effect of the RGB-D camera, a periodically repeated expected motion trajectory was used, while the MARG sensor used a random expected trajectory. The experimental results are shown in Figures8 and 9.Figure 8
Results of depth image enhancement.
(a)(b)(c)Figure 9
Results of attitude calculation enhancement.Comparing Figures8(a) and 8(b), it can be seen that the point cloud after enhancement is closer to the real environment and has a less invalid blind area caused by the narrow environment. Figure 8(c) further demonstrates the impact of depth image enhancement on the accuracy of position estimation. It can be seen from the figure that the position trajectory obtained by using the original image information fluctuates more drastically and deviates from the reference trajectory to a greater degree. But the position trajectory obtained by using the enhanced image information is closer to the reference trajectory, and the overall fluctuation is smaller, which is more conducive to the fusion pose estimation with other sensors.Comparing the desired and enhanced quaternion attitude in Figure9, the overall calculation error of the QEKF algorithm proposed in this paper is within the expected range, which means acceleration information can compensate for the integral drift of angular velocity and geomagnetic field information can calibrate angle measurement errors. Through the enhancement of attitude calculation, it can provide more reliable attitude data for multisensor fusion.


## 5.3. Simulation of the Fusion Part

Because the statistical tests are very sensitive when the sample size is small, both simulation tests and actual experiments were carried out to complete the verification of the multisensor fusion pose estimation method. Simulation tests were mainly used to test the stability of the method under extreme conditions, and actual experiments were mainly used to test the fusion effect of the method in real environments. It was easy to adjust the sensor parameters in simulation environments, so a large number of samples can be generated to test the performance of the method. Although the experimental samples that were collected in the real environments were limited, these samples were more targeted and reliable after the screening of the experimental methods by the simulation tests. At the same time, to present the experimental results as objectively as possible, the root mean square (RMS) errors of multiple experimental results were used to evaluate the performance differences of different methods.Since real sensors can only be simulated numerically under simulation conditions, the actual state and sensors were simplified during the simulation. Different from the real system, an observation state and two sensors were set up in the simulation environment. The observation noise of sensor 1 was smaller than that of sensor 2, and sensor 1 only produced a stable zero offset, while the zero offset of sensor 2 increased with time. Figure10 and Table 1 show the effect of the synchronization update method on the fusion result.Figure 10
The impact of synchronization.
(a)
Single sensor 1(b)
Single sensor 2(c)
Sensor 1 and 2 update frequency:f1>f2(d)
Sensor 1 and 2 update frequency:f1<f2(e)
Sensor 1 and 2 synchronizedTable 1
RMS errors of different update methods.
Update methodSensor1SingleSensor2SingleSensors 1 and 2f1>f2Sensors 1 and 2f1>f2Sensors 1 and 2SynchronizedRMS error1.021.160.8240.9620.413Since the noise of sensor 2 is greater than that of sensor 1, the RMS error of sensor 2 is greater than that of sensor 1 when the state of a single sensor is estimated. In the simulation, the data update frequencyf1 of sensor 1 was different from the data update frequency f2 of sensor 2. If the data synchronization of the sensor was not considered, the state update was performed when the data of any sensor inputted. From the experimental results shown in Figures 10(c) and 10(d), the RMS error when f1>f2 is less than that when f1<f2, indicating that the fusion result is affected by the sensor update frequency, and it is biased towards high-sampling-rate data sources. If the high-sampling-rate data is noisy, the fusion accuracy cannot be significantly improved. The above defects are improved in the modified synchronization update method. It can be seen from the experimental results shown in Figure 10(e) that the RMS error of the fusion data is significantly reduced, and there will be no sawtooth fluctuations in the asynchronous fusion. Therefore, the data synchronization of sensors is an important factor in multirate sensor fusion.In the simulation environment, the process and observation noise conformed to the Gaussian distribution, but the variance of the distribution changed at every moment and satisfied a nonlinear model. The RBF network was designed to predict the variances of the process and observation noise and was trained when the filter was updated. In the simulation, the variances predicted by the network in the test set were compared with the preset variances, and the coefficient of determination (marked asR-square) was calculated. The closer the coefficient of determination is to 1, the higher the accuracy of the network. It can be seen from the results shown in Figure 11 that the network can predict the variances of the noise accurately and has no low coefficients of determination (0.71 for the process noise and 0.84 for the observation noise).Figure 11
Prediction of the noise matrix.Although the coefficient of determination cannot be obtained when the true variance is unknown, the accuracy of the network can be indirectly reflected through the state estimation deviation. In other words, if the state estimation deviation can remain stable even in extreme cases, it means that the network’s prediction of the parameters is accurate. In the simulation, the process noise and the observation noise of the sensor were designed to increase with time. Figure12 and Table 2 show the fusion effects of different fusion algorithms in this case.Figure 12
The fusion effects of different fusion algorithms.
(a)(b)Table 2
RMS errors of different fusion methods.
Fusion methodEKFUKFRBF-AEKFRBF-AUKFRMS error39.4633.4214.987.863It can be seen from the experimental results that at the beginning, the system noise is small and the estimation error of each fusion method is not much different, but as time goes by, the system noise increases. The original EKF and UKF algorithms cannot adaptively adjust the filter parameters, which are easily affected by noise and produce state estimation deviations. Through the RBF network proposed in this paper, the noise covariance matrix can be predicted, so that the EKF and UKF algorithms have adaptive adjustment capabilities. The modified RBF-AEKF and RBF-AUKF methods have smaller RMS errors than the original EKF and UKF methods, so they can be more robust in complex environments.


## 5.4. Experiment of the Fusion Part

Since the environment established in the simulation can only approximately reflect the actual situation, it is still necessary to conduct experiments in the real environment. In actual experiments, the end-effector of the manipulator periodically moved around theZ-axis of the base coordinate system, and the pose of the end-effector was set as the desired value. The position estimations px,py,pz along the X, Y, and Z axes were obtained by fusion of the visual odometer based on the RGB-D camera and the kinematics calculation based on joint encoders. And the quaternion attitude estimations q0,q1,q2,q3 were obtained by fusion of the MARG sensor, visual odometer, and encoders. The experiment compares the original UKF method and the modified adaptive UKF method proposed in this paper, and the results are shown in Figures 13 and 14 and Table 3.Figure 13
Results of position estimation.Figure 14
Results of attitude estimation.Table 3
RMS errors of the original method and modified method.
RMS errorEstimated variablepx(cm)py(cm)pz(cm)q0q1q2q3MARG———0.3190.1010.2600.273VO2.767.022.190.1090.1940.2200.232Encoder0.3031.080.3600.09330.06240.05840.133Fusion (modified)0.1750.3640.1080.008050.005640.006450.0482Fusion (original)0.6631.590.4750.008050.005640.006450.0482It can be seen from Table3 that the RMS error of the pose estimation after multisensor fusion is smaller than that of a single sensor. Compared with that of the original method, the RMS error of the modified method proposed in this paper is significantly reduced in position estimation, but there is little difference in attitude estimation. The main reason is the modified method proposed in this paper is better than the traditional method due to its better adaptive effect when the noise fluctuation is large. In this experiment, the overall fluctuation of the sensor in attitude measurement is smaller than that of position measurement, so the modified method is close to the original method in attitude estimation.It can be seen from Figure13 that the visual odometer is prone to local sudden changes due to environmental interference and the encoder is prone to measurement offset. Compared with the original method, the modified method handles the measurement data of local sudden change better, so the overall RMS error is smaller. And because of the combination of absolute and relative updates, the overall update rate is maintained. In this experiment, the original method was updated only 258 times, while the modified method was updated 479 times, so the real-time response rate of state estimation was improved. It can be seen from Figure 14 that the rotation angle of the desired trajectory on the X and Y axes is small, and the rotation angle on the Z-axis is large. When the rotation angle is small, it can be seen that the sensor has obvious measurement drift, when the rotation angle is large, it is easy to produce local sudden changes. The modified method in this paper can compensate the measurement drift of each sensor, reduce the influence of local sudden changes on the state estimation, and thus can better follow the desired trajectory.


## 6. Conclusions

In order to improve the accuracy of the end-effector pose estimation, this paper proposes a modified method based on multisensor fusion, including preenhancement and adaptive fusion. The following conclusions can be drawn from the experimental results. Firstly, the depth image obtained by the RGB-D camera and the attitude obtained by the MARG sensor is enhanced, which can make the multisensor fusion process in a complex environment more stable. Secondly, the synchronization of sensors has a great influence on the results of multisensor fusion. The control variable is relatively updated when asynchronous, and the state variable is absolutely updated when synchronous, which can balance the update frequency and fusion accuracy. Thirdly, predicting the noise matrix through a neural network can make the algorithm adjust the fusion ratio of each sensor adaptively according to changes in the environment, which helps improve robustness. Finally, in the simulation environment, the RMS errors of the original UKF and EKF method are 4.25 times and 5.02 times that of the modified RBF-UKF method, respectively, and in the actual environment, the average pose estimation errors of the single sensor and the original UKF method are 8.57 times and 3.91 times that of the modified RBF-UKF method, and the update frequency of the original method is only 54% of the modified method. In summary, the method proposed in this paper can improve the accuracy of end-effector pose estimation and has high application value and practical engineering value.



---

*Source: 5550850-2021-04-16.xml*