An Adaptive UWB / MEMS-IMU Complementary Kalman Filter for Indoor Location in NLOS Environment

: High precision positioning of UWB (ultra-wideband) in NLOS (non-line-of-sight) environment is one of the hot issues in the direction of indoor positioning. In this paper, a method of using a complementary Kalman ﬁlter (CKF) to fuse and ﬁlter UWB and IMU (inertial measurement unit) data and track the errors of variables such as position, speed, and direction is presented. Based on the uncertainty of magnetometer and acceleration, the noise covariance matrix of magnetometer and accelerometer is calculated dynamically, and then the weight of magnetometer data is set adaptively to correct the directional error of gyroscope. Based on the uncertainty of UWB distance observations, the covariance matrix of UWB measurement noise is calculated dynamically, and then the weight of UWB data observations is set adaptively to correct the position error. The position, velocity and direction errors are corrected by the fusion of UWB and IMU. The experimental results show that the algorithm can reduce the gyroscope deviation with magnetic noise and motion noise, so that the orientation estimates can be improved, as well as the positioning accuracy can be increased with UWB ranging noise.


Introduction
Ultra-wideband (UWB) positioning technology in indoor positioning can achieve decimeter positioning accuracy, and has been widely used in hospital patient tracking, UAV (unmanned aerial vehicle) positioning, supermarket shopping and so on [1][2][3]. In a line-of-sight (LOS) environment, UWB has high ranging accuracy. The literature [4] shows that the accuracy is about 320 ± 30 mm in an indoor office room. However, in some special scenarios, such as warehousing robot location, emergency rescue, etc., the indoor environment is more complex, and the UWB signal may be blocked by people, goods, or other obstacles, resulting in the signal multi-path effect, the intensity attenuation, even signal loss and so on, which leads to the sharply decline of UWB positioning accuracy [5][6][7]. How to obtain high positioning accuracy in the case of non-line-of-sight (NLOS) environment is a hot issue for researchers [8][9][10]. The fusion of inertial measurement unit (IMU) and UWB is a trend to realize indoor high precision and real-time positioning. Through the IMU integral data, the observations of speed, direction, and position can be obtained. To a certain extent, it can not only eliminate the CKF has three advantages [20,30,31]: first, there is no need to define a motion model, that is, the algorithm can be used for both vehicle and pedestrian applications. Secondly, the error of the state is stored rather than the state itself. When the system is linearized, it can be approximated by a smaller quantity, and a relatively more accurate result can be obtained. Finally, after each correction of the state value, it can be assumed that the current state is correct, that is, the error of the state value is reset to zero. In this way, in each prediction calculation of CKF, only the calculation of system covariance P is needed. In the update stage, the measurement matrix H can be simplified because the observation equation has no prior error.
The remainder of the paper is organized as follows: In Section 2, the fusion algorithm of UWB and IMU based on CKF is discussed, and the motion model and observation model of the algorithm are given. The uncertainty setting method of different observation values in the observation model is analyzed in detail. Subsequently, several experiments are analyzed in Section 3, and Section 4 concludes the paper.

State Model
Following the method in [32,33], the state of CKF algorithm is defined as follows: where is the position error, is the velocity error, is the directional error, is the bias of the gyroscope, and is the bias of acceleration. is a 15-dimensional vector. The first three navigation state vectors are under the navigation coordinate system (n-frame), and the last two bias vectors are under body frame (b-frame). Before integrating the IMU data, the bias vector needs to be removed.
The differential equation of system dynamic model under continuous time is defined as follows: = + (1) is the state transition matrix and is the system noise. In order to facilitate the processing of EKF algorithm, the formula (1) is discretized into = ∅ + (2) ∅ = (3) is the transition matrix at k time. In order to determine , the transformation formula of state X must be derived. CKF has three advantages [20,30,31]: first, there is no need to define a motion model, that is, the algorithm can be used for both vehicle and pedestrian applications. Secondly, the error of the state is stored rather than the state itself. When the system is linearized, it can be approximated by a smaller quantity, and a relatively more accurate result can be obtained. Finally, after each correction of the state value, it can be assumed that the current state is correct, that is, the error of the state value is reset to zero. In this way, in each prediction calculation of CKF, only the calculation of system covariance P is needed. In the update stage, the measurement matrix H can be simplified because the observation equation has no prior error.
The remainder of the paper is organized as follows: In Section 2, the fusion algorithm of UWB and IMU based on CKF is discussed, and the motion model and observation model of the algorithm are given. The uncertainty setting method of different observation values in the observation model is analyzed in detail. Subsequently, several experiments are analyzed in Section 3, and Section 4 concludes the paper.

State Model
Following the method in [32,33], the state of CKF algorithm is defined as follows: X = [ δp n δv n ε b g b a ], where δp n is the position error, δv n is the velocity error, ε is the directional error, b g is the bias of the gyroscope, and b a is the bias of acceleration. X is a 15-dimensional vector. The first three navigation state vectors are under the navigation coordinate system (n-frame), and the last two bias vectors are under body frame (b-frame). Before integrating the IMU data, the bias vector needs to be removed.
The differential equation of system dynamic model under continuous time is defined as follows: . F is the state transition matrix and W is the system noise. In order to facilitate the processing of EKF algorithm, the formula (1) is discretized into F k is the transition matrix at k time. In order to determine F k , the transformation formula of state X must be derived.

The Equation of Acceleration and Gyroscope Bias
The measuring equations of acceleration and gyroscope are as follows: f b and ω b are the measured values of acceleration and gyroscope, f b and ω b are the true values of acceleration and gyroscope, n a and n g are the measurement noise of the acceleration and the gyroscope, respectively, which satisfy the Gaussian distribution, and the covariance are defined as N a ,N g , respectively. b a and b g are drift bias of the acceleration and the gyroscope, which are not a static value, but a time-dependent first-order Markov process, defined as follows: µ a and µ g are offset noise of the acceleration and the gyroscope, respectively, and satisfies Gaussian distribution, and the covariance are defined as U a , U g , respectively.

The Directional Error Equation
The directional error ε is caused by the gyro error and is defined as follows: C n b represents the conversion from b-frame to n-frame, while δω b is the measurement error of the gyroscope, which is caused by drift bias and noise:

The Error Equation of Velocity and Position
The velocity error is caused by acceleration error and direction error, and is defined as: [ f n ×]ε represents the effect of directional error on acceleration, δ f b is the measurement error of acceleration, which is caused by drift bias and noise and is defined as: The position error is defined as: δ From the error equation, the state transition matrix F of the system is deduced as follows: The system noise W is: The covariance matrix of system noise W is: Define the noise transformation matrix as: The noise covariance matrix is discretized to: Thus far, ∅ k , Q k in the prediction equation of CFK algorithm was defined, while Z k , H k , R k and other matrices were defined in Section 2.2.

UWB Observations
There are n base stations with known coordinates (four used in the experiment), The observation function is defined as: wherep ins,k = (x k ,ŷ k ,ẑ k ) represents the position coordinates calculated by INS. δp n k = (δx k , δŷ k , δẑ k ) represents a priori position error calculated by the state transition equation. . represents the Euclidean distance. The formula (18) represents the difference between the ranging value after considering the Observations are defined as: r k is the ranging data of UWB. The formula (20) represents the distance difference between the ranging data of the UWB and the distance to the base station calculated from the INS position.
The standard observation equation is defined as follows: For UWB data, the accuracy of UWB sensor observations may be affected by ambient temperature, power supply stability, fixed obstacles, and even people or objects moving in the positioning scene. Therefore, it is necessary to estimate the confidence of UWB observations. A correlation method based on UWB positioning residual is designed.
In order to simplify the description, there are three base stations, and in the plane location, it is assumed that the UWB base station with the coding number i is recorded as Beacon i , and its seat is marked (x i , y i ). The UWB tag used for positioning is marked tag, and its seat is marked (x, y). The real distance between tag and Beacon i is recorded as r i , and the measured value of this distance is recorded as r i . As Figure 2 shows, ideally, r i = r i , where the only point can be determined by the intersection of three circles, is the location of the tag under the current observation data. In order to obtain the solution of this point, an error function is defined, and the coordinates of tag are obtained by minimizing the error function. One possible error function is: where, abs() represents an absolute function. The coordinates (x , y ) can be obtained by minimizing E(x, y), that is: Ideally, the minimum value of E(x, y) is 0. However, in practical application, the measured value often has a certain deviation. Suppose r i still represents the real distance, and the deviation of r i is recorded as ∆r i . At this point, the measured value of the corresponding Beacon i is r i = r i + ∆r i . Figure 3 shows the area where the circles represented by the three beacon measurements around the tag intersect. In the case of errors in the measurements, the three circles do not intersect at one point, but intersect in pairs. In this case, the estimate of the tag coordinates (x , y ) is still obtained by minimizing E(x, y). Generally speaking, when the error of r i is large, the value of ∆r i is also larger.
The uncertainty for defining UWB data is as follows: C ui ∈ [0, +∞], i ∈ [1, n]. The value of α u is a constant, which is an estimate of the overall residual value and indicates the degree of distrust in the original UWB data. When the value of α u or ∆r i is higher, the more unreliable the data is. C u_max is to limit the value of C ui in a certain range, otherwise it is easy to appear singular matrix. The measurement noise matrix R k of UWB is defined as follows: σ 2 si represents the covariance of the ranging value from the No.i UWB base station to the tag. If there are less than three values in the ranging data due to occlusion and other reasons, the effective residual cannot be obtained. In the experiment, if the sum of the two ranging values is greater than the distance between the two base stations, the residual is 0, and C u = α u . The range values that are obscured and cannot be obtained should be removed from the corresponding rows in R k and Z k .
intersect. In the case of errors in the measurements, the three circles do not intersect at one point, but intersect in pairs. In this case, the estimate of the tag coordinates ( , ) is still obtained by minimizing ( , ). Generally speaking, when the error of is large, the value of Δ is also larger.  The uncertainty for defining UWB data is as follows: The value of α is a constant, which is an estimate of the overall residual value and indicates the degree of distrust in the original UWB data. When the value of α or Δr is higher, the more unreliable the data is. C _ is to limit the value of C in a certain range, otherwise intersect. In the case of errors in the measurements, the three circles do not intersect at one point, but intersect in pairs. In this case, the estimate of the tag coordinates ( , ) is still obtained by minimizing ( , ). Generally speaking, when the error of is large, the value of Δ is also larger.  The uncertainty for defining UWB data is as follows: The value of α is a constant, which is an estimate of the overall residual value and indicates the degree of distrust in the original UWB data. When the value of α or Δr is higher, the more unreliable the data is. C _ is to limit the value of C in a certain range, otherwise

Observations of Accelerometers and Magnetometers
The error observations of accelerometers and magnetometers are defined as follows: f b andm b are offset compensated measurements, respectively.g n and m n are the standard values of gravity acceleration and geomagnetic field in the navigation system.
The observation vector is defined as: The observation equation is defined as follows: The observation model matrix H k is defined as follows: In general, the measurement noise covariance matrix R k is constant, but in fact, R k is closely related to the specific sensor application environment. For gravity data, pitch and roll angle can be estimated. On the other hand, the geomagnetic data can be used to estimate the yaw angle. However in fact, body acceleration cannot be neglected with comparison to gravity, and the geomagnetic data will also have disturbance, which will affect the accuracy of the direction correction. In order to reduce the influence of the direction calculation of the noise data, a method based on dynamic R k matrix is designed.
The magnitude residual (E m m ) of the magnetometer and the magnitude residual (E m g ) of the gravity acceleration are defined as follows: . is taken as the normal value. Both E m m and E m g ∈ [1, +∞], which represent the degree of deviations of the measured valuesm b andf b from the standard values m n and g n on magnitude. The direction residual (E d m ) of the magnetometer and the direction residual (E d g ) of the gravity acceleration are also defined as follows: Both E d m and E d g ∈ [0, +∞]. Define the uncertainty of the magnetometer and acceleration data as follows: Both C m and C g ∈ [0, +∞]. C m_max and C g_max are to limit the values of C m and C g in a certain range, otherwise singular matrices are easy to appear. Both α m and α g are greater than zero, indicating the degree of distrust in the original data of the magnetometer and the acceleration. The greater value of α or E, the less reliable the data is. Finally, the mean of R k is defined as follows: m represent the covariance of the acceleration and magnetometer measurements, respectively.

INS Navigation Equation
Through the integration of gyroscope and acceleration data, the direction, velocity, position, and other data under the navigation system can be obtained. The INS navigation equation is defined as follows: g n is the gravity vector under n-frame, f b is the acceleration vector under b-frame, is the skew-symmetric matrix of the angular velocity, as follows: For acceleration and gyroscope observations at k time, the offset value calculated by the CKF algorithm is first removed), as follows: f b k and ω b k represent the original observations of the acceleration and angular velocity, andf b k and ω b k represent the offset compensated values. The transformation of acceleration from k time to k+1 time is as follows: ⊗ represents the vector cross multiplication and the rotation correction of the acceleration caused by the change of angular velocity.
The transformation of the speed from the k time to the k+1 time is as follows: δv n k represents the error of the velocity at k time, which is cyclically calculated by the CKF algorithm. The change of position from k time to k+1 time is as follows: δp n k represents the error of the position at k time, which is cyclically calculated by the CKF algorithm. The change of posture from k time to k+1 time is as follows: Remote Sens. 2019, 11, 2628 10 of 21 The value of C n b,k+1 represents the result of correcting the orientation error of the pose change matrix C n b,k at k time, and the error ε is calculated cyclically by CKF algorithm. Then, for C n b,k+1 , the rotation compensation, from k time to k + 1 time is carried out to obtain the rotation matrix C n b,k+1 . In order to improve the accuracy of attitude values, the value of C n b,k+1 requires periodic normalization.

Experiment
A test site is set up in the underground garage of a university. As shown in Figure 4, four UWB base stations are placed on four corners of the rectangular area. The IMU device chose X-IMU from X-IO company in the UK. The core chip of UWB tag/base station is DWM1000 chip of Decawave Company. The UWB tag is fixed to the car. At the same time, the IMU is fixed 5 cm below the UWB tag. In order to keep the data of IMU and UWB synchronized, the notebook in the car receives ranging data from IMU and UWB at the same time. In the process of clockwise movement, the car basically maintains a uniform speed.
cannot be obtained, this trajectory is only used as a reference schematic diagram. In addition, in order to increase the location noise of the eight-shaped route, one to two base stations will be blocked randomly, which will cause the ranging data to be abnormal or data loss.  In the experimental analysis, first of all, the positioning results are analyzed. Then the influence on the positioning results is analyzed from different types of observations, such as geomagnetic residual, gravity residual, ranging residual and so on.

Positioning Trajectory Analysis
The location result of path 1 is given in Figure 6. The left figure is the location result given by the optimization algorithm. In general, with the exception of a few points deviating from the overall trajectory, the rest of the anchor points are aggregated on the path of the car motion. The right image shows the trajectory of the fusion of UWB and IMU. It can be seen that the fusion trajectory gives a higher frequency of positioning results, and provides attitude information for each positioning result. Compared with the left image, the track of the right image is smoother. The experiment set two routes, one is a rectangular route with fewer turns, and the other is an eight-shaped route with more turns, such as route 1 and route 2 in Figures 4 and 5. The locations of the four base stations are shown in the rectangular points, and the two red circles in the figure represent the two columns in the underground garage. Because the real trajectory of the car motion cannot be obtained, this trajectory is only used as a reference schematic diagram. In addition, in order to increase the location noise of the eight-shaped route, one to two base stations will be blocked randomly, which will cause the ranging data to be abnormal or data loss.
annot be obtained, this trajectory is only used as a reference schematic diagram. In addition, in order o increase the location noise of the eight-shaped route, one to two base stations will be blocked andomly, which will cause the ranging data to be abnormal or data loss.  In the experimental analysis, first of all, the positioning results are analyzed. Then the influence n the positioning results is analyzed from different types of observations, such as geomagnetic esidual, gravity residual, ranging residual and so on. In the experimental analysis, first of all, the positioning results are analyzed. Then the influence on the positioning results is analyzed from different types of observations, such as geomagnetic residual, gravity residual, ranging residual and so on. The location result of path 1 is given in Figure 6. The left figure is the location result given by the optimization algorithm. In general, with the exception of a few points deviating from the overall trajectory, the rest of the anchor points are aggregated on the path of the car motion. The right image shows the trajectory of the fusion of UWB and IMU. It can be seen that the fusion trajectory gives a higher frequency of positioning results, and provides attitude information for each positioning result. Compared with the left image, the track of the right image is smoother.  In the experimental analysis, first of all, the positioning results are analyzed. Then the influence on the positioning results is analyzed from different types of observations, such as geomagnetic residual, gravity residual, ranging residual and so on.

Positioning Trajectory Analysis
The location result of path 1 is given in Figure 6. The left figure is the location result given by the optimization algorithm. In general, with the exception of a few points deviating from the overall trajectory, the rest of the anchor points are aggregated on the path of the car motion. The right image shows the trajectory of the fusion of UWB and IMU. It can be seen that the fusion trajectory gives a higher frequency of positioning results, and provides attitude information for each positioning result. Compared with the left image, the track of the right image is smoother.  The main state of the fusion algorithm is shown in Figure 7. From top to bottom, they are the bias of the acceleration, the bias of the gyroscope, the change of velocity, and angle in turn. Because the residual of UWB ranging data is relatively small, the integral drift of INS is well constrained, and the bias of velocity and acceleration is effectively modified. As can be seen from subfigure 3 of Figure 7, the speed value is basically about ±1 meter, which is consistent with the actual value. As can be seen from subfigure 4 of Figure 7, through the correction of the opposite direction value of the magnetometer and accelerometer, the bias of the gyroscope can be estimated effectively. In addition, the angle change of the trolley movement on path 1 is relatively regular, and the angle value is close to the actual value. Subfigures 1 and 2 are bias values estimated from position and angle updates, which are subtracted from each integral of the accelerometer and gyroscope. 7, the speed value is basically about ±1 meter, which is consistent with the actual value. As can be seen from subfigure 4 of Figure 7, through the correction of the opposite direction value of the magnetometer and accelerometer, the bias of the gyroscope can be estimated effectively. In addition, the angle change of the trolley movement on path 1 is relatively regular, and the angle value is close to the actual value. Subfigures 1 and 2 are bias values estimated from position and angle updates, which are subtracted from each integral of the accelerometer and gyroscope.

Analysis of Geomagnetism and Gravity Residual
From the above figure, it can be seen that the magnetic field data in the whole path 1 is relatively stable, and the geomagnetism changes regular with the attitude change of the trolley in two consecutive circles. Initially, the car stood still. The geomagnetic data for a period of time (5 to 10 seconds) are collected and converted to n-frame as the environmental standard geomagnetic data ( in the second part of Section 2.2). If there is no geomagnetic interference in the environment, all subsequent geomagnetic data should be equal to the standard geomagnetic data after they are converted to n-frame. On this basis, the transformation matrix from b-frame to n-frame in the whole process of motion can be obtained.

Analysis of Geomagnetism and Gravity Residual
From the above figure, it can be seen that the magnetic field data in the whole path 1 is relatively stable, and the geomagnetism changes regular with the attitude change of the trolley in two consecutive circles. Initially, the car stood still. The geomagnetic data for a period of time (5 to 10 seconds) are collected and converted to n-frame as the environmental standard geomagnetic data (m n in the second part of Section 2.2). If there is no geomagnetic interference in the environment, all subsequent geomagnetic data should be equal to the standard geomagnetic data after they are converted to n-frame. On this basis, the transformation matrix C n b from b-frame to n-frame in the whole process of motion can be obtained.
The bottom figure in Figure 8 shows the difference between the magnetometer data converted to n-frame and the environmental standard geomagnetic data m n , which should ideally be zero. The non-zero value in the graph comes from the influence of two cases: one is that the geomagnetism in the environment is unstable; the other is that the calculation of the transformation matrix C n b itself is not completely accurate. However, due to the relatively small error, it can be seen from Table 1 that the average geomagnetic error is not more than 0.2 g. Therefore, it is theoretically feasible to calculate the direction of geomagnetism in this environment.  Figure 8 shows the difference between the magnetometer data converted to n-frame and the environmental standard geomagnetic data , which should ideally be zero. The non-zero value in the graph comes from the influence of two cases: one is that the geomagnetism in the environment is unstable; the other is that the calculation of the transformation matrix itself is not completely accurate. However, due to the relatively small error, it can be seen from Table 1 that the average geomagnetic error is not more than 0.2 g. Therefore, it is theoretically feasible to calculate the direction of geomagnetism in this environment.

Coordinate axis Average of absolute value Max
Min X 0.0449 0.1188 -0.1180 Y 0.1564 0.2720 -0.0125 Z 0.0287 0.0839 -0.0783 For gravity data, similar to magnetometer data, accelerometer data collected for an initial rest period (5 to 10 seconds) are converted to n-frame as environmental standard gravity data ( in the second part of Section 2.2). If the interference of acceleration in the motion of the car is ignored, all subsequent accelerometer data should be equal to the standard gravity data converted to n-frame. On this basis, the transformation matrix from b-frame to n-frame in the whole process of motion can be obtained.  For gravity data, similar to magnetometer data, accelerometer data collected for an initial rest period (5 to 10 seconds) are converted to n-frame as environmental standard gravity data (g n in the second part of Section 2.2). If the interference of acceleration in the motion of the car is ignored, all subsequent accelerometer data should be equal to the standard gravity data converted to n-frame. On this basis, the transformation matrix C n b from b-frame to n-frame in the whole process of motion can be obtained.
As shown in Figure 9, the above figure shows the acceleration data for path 1, and the following figure shows the difference between the accelerometer data converted to n-frame and the standard gravity data g n , which should ideally be zero. The non-zero value in the graph comes from the influence of two cases: one is that the motion of the trolley is actually non-uniform and moving under the acceleration caused by the periodic pulling force; the other is that the calculation of the transformation matrix C n b itself is not completely accurate. As can be seen from Table 2, the error between X-and Z-axis is relatively small, and the average gravity residual is less than 0.8 G. The average gravity residual of the Y-axis is about 2 G, which explains why the acceleration residual of the Y-axis is increasing in subfigure 1 of Figure 7. Generally speaking, the calculation of gravity constraint direction in this environment is still feasible in theory, but its confidence is lower than that of geomagnetic data. As shown in Figure 9, the above figure shows the acceleration data for path 1, and the following figure shows the difference between the accelerometer data converted to n-frame and the standard gravity data , which should ideally be zero. The non-zero value in the graph comes from the influence of two cases: one is that the motion of the trolley is actually non-uniform and moving under the acceleration caused by the periodic pulling force; the other is that the calculation of the transformation matrix itself is not completely accurate. As can be seen from Table 2, the error between X-and Z-axis is relatively small, and the average gravity residual is less than 0.8 G. The average gravity residual of the Y-axis is about 2 G, which explains why the acceleration residual of the Y-axis is increasing in subfigure 1 of Figure 7. Generally speaking, the calculation of gravity constraint direction in this environment is still feasible in theory, but its confidence is lower than that of geomagnetic data.

Residual Analysis of UWB Data
As can be seen in Figure 10, the UWB ranging data of route 1 is partially obscured by the column, resulting in a state of discontinuity. The corresponding Δ is shown in Figure 11. The columns of different colors represent the ranging residual value of the corresponding base station respectively. As can be seen from Figure 10, the ranging value is relatively unstable at several peaks, because the ranging error of UWB increases with the increase of distance. In Figure 11, the ranging residual increases accordingly. By analyzing the proportion of ranging residual, the weight of covariance of observation values of different base stations is determined dynamically. As can be seen from Table 3, the average ranging residual of the four base stations is less than 0.1 meters, indicating that the overall quality of path 1 is high. The minimum residual value of 0 indicates that the ranging value is obscured at a certain time, and the residual error is not calculated at this time.

Residual Analysis of UWB Data
As can be seen in Figure 10, the UWB ranging data of route 1 is partially obscured by the column, resulting in a state of discontinuity. The corresponding ∆r i is shown in Figure 11. The columns of different colors represent the ranging residual value of the corresponding base station respectively. As can be seen from Figure 10, the ranging value is relatively unstable at several peaks, because the ranging error of UWB increases with the increase of distance. In Figure 11, the ranging residual increases accordingly. By analyzing the proportion of ranging residual, the weight of covariance of observation values of different base stations is determined dynamically. As can be seen from Table 3, the average ranging residual of the four base stations is less than 0.1 meters, indicating that the overall quality of path 1 is high. The minimum residual value of 0 indicates that the ranging value is obscured at a certain time, and the residual error is not calculated at this time.

Location Trajectory Analysis
The location result of path 2 is given in Figure 12, and the left figure is the location result given by the optimization algorithm. Because more positions on path 2 are blocked by stone columns, coupled with the human interference of the testers, there are many of anomalies that deviate from

Location Trajectory Analysis
The location result of path 2 is given in Figure 12, and the left figure is the location result given by the optimization algorithm. Because more positions on path 2 are blocked by stone columns, coupled with the human interference of the testers, there are many of anomalies that deviate from

Location Trajectory Analysis
The location result of path 2 is given in Figure 12, and the left figure is the location result given by the optimization algorithm. Because more positions on path 2 are blocked by stone columns, coupled with the human interference of the testers, there are many of anomalies that deviate from the trajectory. The right image is the trajectory of UWB and IMU fusion, it can be seen that the fusion track gives a higher frequency positioning results, compared with the left image, the right trajectory effectively shielded the abnormal points of the left image, and the trajectory is smoother. the trajectory. The right image is the trajectory of UWB and IMU fusion, it can be seen that the fusion track gives a higher frequency positioning results, compared with the left image, the right trajectory effectively shielded the abnormal points of the left image, and the trajectory is smoother.

Analysis of Geomagnetism and Gravity Residual
The blue data in Figure 13 is the residual modulus value of gravity data ‖ ‖ , and the corresponding red data is the uncertainty of gravity residual data. The blue data in Figure 14 is the residual modulus value of the geomagnetic data ‖ ‖, which corresponds to the uncertainty of the geomagnetic residual data.

Analysis of Geomagnetism and Gravity Residual
The blue data in Figure 13 is the residual modulus value of gravity data δ f n , and the corresponding red data is the uncertainty C g of gravity residual data. The blue data in Figure 14 is the residual modulus value of the geomagnetic data δm n , which corresponds to the uncertainty C m of the geomagnetic residual data. the trajectory. The right image is the trajectory of UWB and IMU fusion, it can be seen that the fusion track gives a higher frequency positioning results, compared with the left image, the right trajectory effectively shielded the abnormal points of the left image, and the trajectory is smoother.

Analysis of Geomagnetism and Gravity Residual
The blue data in Figure 13 is the residual modulus value of gravity data ‖ ‖ , and the corresponding red data is the uncertainty of gravity residual data. The blue data in Figure 14 is the residual modulus value of the geomagnetic data ‖ ‖, which corresponds to the uncertainty of the geomagnetic residual data.   The residual of the magnetic field data and gravity data in the whole path 2 are basically consistent with the uncertainty of setting, which shows that the given uncertainty formula and parameters can well reflect the residual variation of geomagnetism and gravity data. Thus, in the process of filtering, we can intervene the covariance of the observed values in order to obtain a more accurate direction value. The correction and update of the direction value further restrict the bias of the gyroscope. Figure 15 is a comparison of the effects of and on gyroscope data. It can be seen that under the condition of uncertainty intervention, the bias error of the gyroscope is almost twice as small, thus improving the accuracy of angle estimation.

Residual Analysis of UWB Data
As can be seen in Figure 16, the UWB ranging value of route 2 has changed greatly at some time due to the occlusion of stone columns and human interference. The corresponding Δ is shown in Figure 17, and the optimal residual also increases at the ranging jump point. As can be seen from Table 4, although the average ranging residual of the four base stations is about 0.2 m, the highest residual value is more than two meters, which will cause the outlier with large positioning deviation. By dynamically increasing the covariance of ranging outliers, the influence of anomalies on the fusion algorithm can be shielded. The residual of the magnetic field data and gravity data in the whole path 2 are basically consistent with the uncertainty of setting, which shows that the given uncertainty formula and parameters can well reflect the residual variation of geomagnetism and gravity data. Thus, in the process of filtering, we can intervene the covariance of the observed values in order to obtain a more accurate direction value. The correction and update of the direction value further restrict the bias of the gyroscope. Figure 15 is a comparison of the effects of C g and C m on gyroscope data. It can be seen that under the condition of uncertainty intervention, the bias error of the gyroscope is almost twice as small, thus improving the accuracy of angle estimation. The residual of the magnetic field data and gravity data in the whole path 2 are basically consistent with the uncertainty of setting, which shows that the given uncertainty formula and parameters can well reflect the residual variation of geomagnetism and gravity data. Thus, in the process of filtering, we can intervene the covariance of the observed values in order to obtain a more accurate direction value. The correction and update of the direction value further restrict the bias of the gyroscope. Figure 15 is a comparison of the effects of and on gyroscope data. It can be seen that under the condition of uncertainty intervention, the bias error of the gyroscope is almost twice as small, thus improving the accuracy of angle estimation.

Residual Analysis of UWB Data
As can be seen in Figure 16, the UWB ranging value of route 2 has changed greatly at some time due to the occlusion of stone columns and human interference. The corresponding Δ is shown in Figure 17, and the optimal residual also increases at the ranging jump point. As can be seen from Table 4, although the average ranging residual of the four base stations is about 0.2 m, the highest residual value is more than two meters, which will cause the outlier with large positioning deviation. By dynamically increasing the covariance of ranging outliers, the influence of anomalies on the fusion algorithm can be shielded.

Residual Analysis of UWB Data
As can be seen in Figure 16, the UWB ranging value of route 2 has changed greatly at some time due to the occlusion of stone columns and human interference. The corresponding ∆r i is shown in Figure 17, and the optimal residual also increases at the ranging jump point. As can be seen from Table 4, although the average ranging residual of the four base stations is about 0.2 m, the highest residual value is more than two meters, which will cause the outlier with large positioning deviation.
By dynamically increasing the covariance of ranging outliers, the influence of anomalies on the fusion algorithm can be shielded. As can be seen from Figure 18, the sample points with large residual errors in Figure 17 also have greater uncertainty values in Figure 18. That is, Figure 18 exponentially magnifies the residual in Figure 17 to increase the sensitivity to the residual data.   As can be seen from Figure 18, the sample points with large residual errors in Figure 17 also have greater uncertainty values in Figure 18. That is, Figure 18 exponentially magnifies the residual in Figure 17 to increase the sensitivity to the residual data.   As can be seen from Figure 18, the sample points with large residual errors in Figure 17 also have greater uncertainty values in Figure 18. That is, Figure 18 exponentially magnifies the residual in Figure 17 to increase the sensitivity to the residual data.  0.2304 9.5135 0 Generally speaking, the fusion algorithm can dynamically give the corresponding confidence according to the gravity residual, geomagnetism residual and ranging residual in route 1 and 2, so as to evaluate the observed values effectively and improve the accuracy and robustness of the fusion algorithm.

Conclusions
In this paper, a fusion location method of UWB and IMU based on adaptive CKF is presented. The algorithm optimizes the location results by adaptively setting the weights of UWB data, magnetic field data, and acceleration data and so on. Firstly，through the analysis of UWB ranging data, it can be seen that in non-line-of-sight environment, the accuracy of UWB ranging will be affected. However, by dynamically increasing the covariance of ranging outliers, the influence of anomalies on the fusion algorithm can be shielded. Secondly, with the analysis of geomagnetic and gravity data in the experimental scene, the residual of geomagnetic and gravity data is small, which can meet the needs of heading calculation. However, the credibility of geomagnetic data is better than that of gravity data. Then, the fusion of UWB and IMU not only can eliminate the multi-path and non-lineof-sight effect caused by occlusion of UWB signal, but also restricts the drift of INS, corrects the deviation of velocity and acceleration, and realizes the complementary advantages. However, there are two good prerequisites in the experiment. That is, the motion of the car is relatively stable and the experimental geomagnetic data are relatively stable. Thus, in the face of more complex motion patterns such as rapid acceleration and stop, as well as more and stronger geomagnetic interference environment, better motion models and algorithms need to be studied and further tested.  Generally speaking, the fusion algorithm can dynamically give the corresponding confidence according to the gravity residual, geomagnetism residual and ranging residual in route 1 and 2, so as to evaluate the observed values effectively and improve the accuracy and robustness of the fusion algorithm.

Conclusions
In this paper, a fusion location method of UWB and IMU based on adaptive CKF is presented. The algorithm optimizes the location results by adaptively setting the weights of UWB data, magnetic field data, and acceleration data and so on. Firstly, through the analysis of UWB ranging data, it can be seen that in non-line-of-sight environment, the accuracy of UWB ranging will be affected. However, by dynamically increasing the covariance of ranging outliers, the influence of anomalies on the fusion algorithm can be shielded. Secondly, with the analysis of geomagnetic and gravity data in the experimental scene, the residual of geomagnetic and gravity data is small, which can meet the needs of heading calculation. However, the credibility of geomagnetic data is better than that of gravity data. Then, the fusion of UWB and IMU not only can eliminate the multi-path and non-line-of-sight effect caused by occlusion of UWB signal, but also restricts the drift of INS, corrects the deviation of velocity and acceleration, and realizes the complementary advantages. However, there are two good prerequisites in the experiment. That is, the motion of the car is relatively stable and the experimental geomagnetic data are relatively stable. Thus, in the face of more complex motion patterns such as rapid acceleration and stop, as well as more and stronger geomagnetic interference environment, better motion models and algorithms need to be studied and further tested.