Adaptive Unscented Kalman Filter for Target Tacking with Time-Varying Noise Covariance Based on Multi-Sensor Information Fusion

In this paper, an innovative optimal information fusion methodology based on adaptive and robust unscented Kalman filter (UKF) for multi-sensor nonlinear stochastic systems is proposed. Based on the linear minimum variance criterion, this multi-sensor information fusion method has a two-layer architecture: at the first layer, a new adaptive UKF scheme for the time-varying noise covariance is developed and serves as a local filter to improve the adaptability together with the estimated measurement noise covariance by applying the redundant measurement noise covariance estimation, which is isolated from the state estimation; the second layer is the fusion structure to calculate the optimal matrix weights and gives the final optimal state estimations. Based on the hypothesis testing theory with the Mahalanobis distance, the new adaptive UKF scheme utilizes both the innovation and the residual sequences to adapt the process noise covariance timely. The results of the target tracking simulations indicate that the proposed method is effective under the condition of time-varying process-error and measurement noise covariance.


Introduction
There have been increasing demands for developing robust, adaptive, and accurate multi-sensor information filters (MSIF), which have been widely applied to many fields such as navigation systems, modern industries, military threat detection, target tracking, and remote sensing [1,2]. Especially in recent years, there have been many state estimation problems in which the processes were often non-linear and uncertain for tracking and navigation, for example [3]. Hence, the theory has been researched and broadly applied to many realistic systems. A method with both adaptability and robustness cannot be realized in real time for a single sensor/observer system [4]. By using a multi-sensor structure, an information fusion algorithm can obtain much more accurate estimations than a single one [5,6]. However, to the best of our knowledge, few methods focus on both adaptability and robustness in MSIF. From a system point of view, there are mainly two different methods to process the data from a multi-sensor [7]. The first method is the centralized filter where all raw sensor data are fed to a central site for processing [1]. The second one is the decentralized or distributed filter where the process is divided between some local filter concurrently to obtain individual raw data-based estimates and one master/center filter to fuse those local estimates to provide a much accurate global optimal estimate [6,8].
The centralized filter is also called measurement fusion, because all observations are directly fused to obtain a final estimate. The main advantage is high accuracy due to the the maximum likelihood criterion for the proper choice of filtering weight. They argued that the method is efficient by adapting the matrices R and/or Q [19]. The basic idea of [19] is to adapt the R by innovation sequences and the Q by residual ones. However, the innovation and residual sequences obtained from the filter are not independent [15].
In [20], based on an adaptively robust EKF, Yuan et al. proposed a PDR/UWB (Pedestrian Dead Reckoning and Ultra-Wide Band) integrated navigation algorithm. To obtain the adaptability, the algorithm takes the positioning scene and the heading as constraints. The robustness of the algorithm is achieved by adopting the idea similar to [6]. However, in many other applications such as tracking and remote sensing, the constraints are not always available or difficult to implement because of the great complexity and the high computational cost [21]. In [22], to handle measurement outliers, the robust estimation reduces the weight of the observation; to handle the model error, an adaptive factor is introduced to balance the adverse effect. This algorithm has inspired many scholars. In [23], two novel quantitative nonlinear observability measures are proposed to get an optimal filter design. However, both [22] and [23] inevitably use the state-dependent calculations to adjust the measured values or the weight matrix, so the same problem as in [20] mentioned above exists. In [24], an adaptive filtering method was proposed. The authors used the residual error to construct a low-pass filter together with the process noise covariance. They argued that the high process noise could be effectively suppressed. However, if measurement noise covariance is not estimated accurately in a timely manner, the residual error will be inaccurate; then, in the next iteration, the residual error may not be adjusted only by the algorithm. So, it is necessary to estimate the measurement noise covariance as accurately as possible. Moreover, it is better that the estimation of the measurement noise covariance is independent of the process noise or estimation of the state [15].
If the estimation of the noise covariance matrix R could be accurately made, it would be a solid foundation to give the matrix Q a relatively accurate estimation. Zhang et al. developed a measurement-based adaptive Kalman filtering algorithm (MAKF) that overcame the instability drawback of improved Sage-Husa adaptive filter for the integrated navigation system [25,26]. Realizing the limitation of MAKF is that the following assumption could not always be held-one of the measurement noise covariances is relatively smaller than the other-the group of Zhang afterwards developed an improved method named redundant measurement noise covariance estimation (RMNCE) [27], which can estimate the noise variance of the measurement and is not affected by the process state estimation error. So, in this work, we utilize the RMNCE, which can deal with the unknown noise covariance in real time, to calculate the measurement noise covarianceR k for each sensor.
In our proposed method, the matrixR k of each sensor is also calculated through innovation sequences, denoted as R in , as mentioned in [28]. We denote the ratio between R in andR k as the indicator to reflect whether there would be non-ignorable process error or not. Furthermore, if the indicator or the hypothesis test theory based on chi-square suggests the existence of process error, the trigger for adaptation is on. To the best of our knowledge, the combination of the two criteria is firstly introduced by this paper.
For a given MSIF problem, without loss of generality, the statistical properties of measurement noise are not reliable, though they could be obtained in advance. So, we adopt the RMNCE method to estimate the variance of measurement noise in multi-sensor system. Additionally, taking all the above adaptive Q estimation algorithms into consideration, a new Q estimation algorithm based on both innovation and residual sequences is given, which is inspired by [16]. Finally, the decentralized architecture is used to fuse the estimations from local sensors.
The contributions of the paper are twofold. Firstly, an efficient algorithm is proposed for the MSIF to detect the process-error based on the indicator, which is combined by hypothesis test theory with the Mahalanobis distance of innovation sequences and the RMNCE. Then, an innovative Q estimation algorithm is proposed by using both innovation and residual sequences. Secondly, to the best of our knowledge, the RMNCE algorithm together with a weighted factor is first introduced into MSIF. To begin with, the covariance of measurement noise obtained by RMNCE is not only used as the measurement noise covariance estimation of each sensor but also as the element to calculate the weighted factor. Then, a novel method is proposed to simplify the calculation of the weighted factor as an alternative to optimal matrix weights in [1]. Moreover, an indicator is also proposed based on the RMNCE to detect whether the process-error exists or not. At last, the simulation results demonstrate that the proposed scheme can increase the tracking precision with both adaptivity and robustness.
The remainder of this paper is organized as follows: in the next section, we describe the standard UKF, the RMNCE, an innovative adaptive UKF (AUKF) proposed by this paper and the decentralized MSIF. Section 3 introduces the adaptive multi-sensor information fusion algorithm (RAUKF-MSIF). Section 4 provides the simulation and discussions. Section 5 finally draws the conclusions.

The Decentralized MSIF and the Proposed Adaptive UKF
Considering a discrete time nonlinear stochastic system with l sensors, the process and measurement models can be described as where X k ∈ R n×1 denotes the state vector, Z i k ∈ R m i ×1 is the measurement collected by sensor i at sampling time instant k, W k−1 and V i k are uncorrelated zero-mean Gaussian white noise with compatible dimension [4,16,25,26,28,29], f (•) and h(•) are the known time-varying nonlinear state transition and measurement function, respectively. Γ k−1 is the system noise-driven matrix with compatible dimension.
The statistical properties assumed about noise processes can be summarized as [4] where δ k,j denotes the Kronecker delta function.
The following assumptions are also made as initial value where the initial state X 0 is independent of W k−1 and V i k , P 0 is the initial estimation error covariance matrix.
In the following Section 2.1, the procedure of the standard UKF is briefly reviewed. In Section 2.2, the RMNCE is introduced. To improve the adaptivity of the estimation, the process-error should be detected timely, so an innovative adaptive UKF is proposed in Section 2.3. In Section 2.4, the structure of MSIF adopted by our work is given.

The Standard UKF
The UKF uses the fact that it should be easier to estimate a nonlinear distribution than to give an approximation of a nonlinear system [29]. In the standard UKF, to generate the sigma points to undergo the nonlinear transformation and calculate the first two moments of the transformed set, the UT, a deterministic sampling technique, is implemented. For the sake of simplicity, only one sensor of (1) is taken into consideration, and the general procedures are as follows: Sensors 2021, 21, 5808 5 of 22 Step 1: Initialization.
whereX 0 is the initial state, P 0 is the initial estimation error covariance.
where n denotes the dimension of the state; λ = α 2 (n + κ) − n is the composite scaling parameter that is used for fine tuning, α is set to 0 ≤ α ≤ 1 and a good default setting on κ is κ = 0 [15].
whereX k/k−1 is the predicted state mean; and P XX is the predicted state covariance. ω where β ≥ 0 is used to incorporate the higher order information of the distribution, according to [15], for Gaussian distribution β = 2 [29].
Step 6: State estimation and error covariance matrix update.
Step 7: Iterate from steps 2 to 6 until all samples are completed. Under the condition of time-varying process-error and measurement noise covariance, we can infer that if Q k−1 in (6) and R k in (9) could not be estimated timely, then inaccurate estimates would be made because, in (10),X k is influenced by K k which is related to R k andX k/k−1 .

Adaptive R Estimation
As proved by Li et al. in [27], for simplicity, assuming Z 1 (k) and Z 2 (k) are independent redundant measurements of a signal Z(k) from two sensors, they can be modeled as where Z T (k) denotes the true value of Z(k), S 1 (k) and S 2 (k) are steady items of the measurement errors, V 1 (k) and V 2 (k) are uncorrelated zero-mean random white noise. The first-order-self-difference (FOSD) ∆Z 1 , ∆Z 2 and the second-order-mutual-difference (SOMD) ∆Z 12 are defined as Under the condition that the sampling interval is short enough, The covariance of the random noise for measurement Z 1 (k) and Z 2 (k) can be estimated as The mathematical expectations in (13) are calculated as follows, because the statistical characteristics are stable over a relatively short period.
where E[∆Z i (k)], i = 1, 2 and E[∆Z 12 (k)] can be calculated by where M is the sliding window width which can be empirically set to 30~60 [25]. In practice, in order to capture the transient behavior of the variation of R k in time while considering the smoothness, a fading memory calculation is implemented as [30] where b is the fading factor and in our work is set to 0.980, R k represents the direct output in (13) at time k,R k is the final covariance matrix.

The Proposed Adaptive UKF Based on the Mahalanobis Distance of the Innovation Vector
As emphasized by many researchers [15,17,31], the parameters of the system model and the distribution of the measurement noise and the process noise could not be maintained as constants in practice all the time. In order to prevent the estimation from deteriorating or even diverging, caused by the model error, it is vital to determine and correct the mismatch between the real process error and the parameters' preset. Herein, in our proposed method, based on hypothesis testing theory [14], the Mahalanobis distance of innovation vector is employed as a criterion to identify whether the system modeling error exists or not [17].
For the sake of simplicity, only one sensor of (1) is taken into consideration, thus, the superscript i is dropped. This paper develops a new method to improve the adaptability of the classical UKF against process model error based on the Mahalanobis distance. However, if the measurement noise matrix R k is also contaminated, then one single sensor cannot cope with this case. Thanks to the RMNCE, in MSIF, the R k can be estimated and is immune to the state estimation.
Define the innovation sequence ε k according to [15] as Then ε k should be zero-mean ; the square of the Mahalanobis distance of the innovation should be χ 2 distributed [4], where m is the degree of freedom. Based on the hypothesis testing theory, let α be the given significance level, and then we have where Pr(•) stands for the probability of a random event, χ 2 m,α denotes the α-quantile of the distribution χ 2 m . If (20) does not hold, it can be deduced with high probability (1 − α) that there exists process-error in the system (1), assuming that the observations are within reasonable bounds. Because there is more than one sensor in the system, (20) should be calculated for each sensor.
Different from [4,17], instead of a single parameter that acts on P k/k−1 We present a new algorithm that takes the advantage of sequences of both innovation and residual to correct the matrix Q directly by a diagonal matrix, and the procedure is as follows. According to [15], the residual η k can be defined as According to [18], based on the principle of orthogonality [12], the residual sequences are uncorrelated from the measurements. So, the following equations should hold: If P k is replaced by its definition in (10), then The traces of both sides should be equal: where H k = ∂h ∂X X=X k , and P * k/k−1 is the predicted covariance without the additive process noise. If (23) does not hold, there will be some change in matrix Q. So, one scalar called the adaptive fading factor is generated based on the Equation (23) [16]: Instead of one scalar being calculated to tune the matrix Q, much higher accuracy could be obtained by using a matrix Λ = diag{λ 1 , λ 2 , . . . , λ m }, in which m is the dimension of the Q.
So, an innovative adaptive method is proposed as To obtain (28), firstly, the relationship of are given; then the matrix Λ k is solved but not in the form of the trace of a matrix. The covariance of innovation is: The covariance of difference sequences between innovation and residual is A deep look should be taken at ε k and η k before calculating (30) where X k/k−1 and X k represent the prediction error and the estimation error, respectively. The first two items in the right side of (30) can be derived based on (31) as Because the sample frequency is high, the Jacobian matrices Then, based on (32) and (33), namely, According to (29), (25) can be rewritten without trace calculation based on (32), (33), and Hence, (28) is obtained by rearranging (35). Normally, the covariance equations stacked above are calculated as

Remark 1.
In Equation (28), the inverse matrices of (H k Γ k−1 ) −1 and Q −1 k−1 are required to be calculated. In practice, Q k−1 is normally a positive definite diagonal matrix, but H k Γ k−1 may be not an invertible matrix, or not even a square matrix. Considering that Λ k is a diagonal matrix, we can obtain the elements by solving the following equation: Remark 2. In Equation (37), H k and H T k are the Jacobian matrices. They normally can be considered as the slope atX k . During the course of numerical calculation, H k and H T k may be very small, so one or more elements of Λ k will be large enough to cause P k to be a negative definite in the next step. To solve this problem, a threshold value is used to limit the element of Λ k . Meanwhile, if λ i is not positive in numerical applications, it is always reset to the absolute value of its estimate [28].
Remark 3. In (28) or (37), the result of the first three items in the right side of equations may be much greater than the last item, because the assumption that the innovation and residual sequences are orthogonal is statistically ideal. In other words, for a limited set of the sample, (36) is commonly biased [28]. However, in many real-time or online systems, the sliding window could not be set very large to meet the assumption. So, the tradeoff should be made between the numerical stability and real-time performance.

The Decentralized MSIF
Due to the expensive computational cost for high-dimension matrices and low stability when some measurements are abnormal [10], i.e., one or more sensors provide information with large noises, or just noises, the centralized architecture is not adopted widely. This work adopts the decentralized MSIF structure similar to [1], which is derived from [32]. For simplicity, the time instant k in X k is dropped in this section [1].
. . , l be unbiased estimators of X in (1) and the estimation errors be Assume that X i and X j are correlated, the covariance and cross covariance are P ii and P ij , respectively. The optimal fusion estimatorX op with matrix weights can be described aŝ where the optimal matrix weights A i , i = 1, 2, . . . , l are to be determined.
The globally optimal information fusion Kalman filterX op of the state X, based on the principle of linear minimum variance, will satisfy the following conditions [1,17]: Condition 1.X op must be the unbiased estimation of X, namely, E X op = E[X]. Condition 2.X op makes tr{P} minimum, in which P is the error covariance matrix ofX op .
For simplicity, we denote A = [A 1 , A 2 , . . . , A l ] T . Our aim is to find A to construct the unbiased estimatorX where A i , i = 1, 2, . . . , l are arbitrary matrices. IfX op is an unbiased estimator for X, the following condition should be fulfilled, Taking the expectation of both sides of (38) yields, So (43) is obtained Based on (38) and (42), we get the fusion estimation error The error variance matrix of the fusion estimator is where Σ = P ij i, j = 1, 2, . . . , l is a symmetric positive definite matrix. Now the problem is converted to a classic one: under the constraint of (42), to solve the minimum of tr{P} = tr A T ΣA by applying the Lagrange multiplier where Λ ∈ R n×n is the Lagrange multiplier, E = [I n , I n , . . . , I n ] T ∈ R nN×n ; and I n represents an n-dimensional identity matrix. By setting ∂F/ ∂A| A=A = 0, and noting that Σ T = Σ, we have By substituting (42) into (46) Because Σ is a symmetric positive definite matrix, E T Σ −1 E is nonsingular. Based on the matrix theory, (47) can be solved From (48), the matrix Σ should be calculated to obtain the global optimal state estimation. The diagonal elements P ii , i = 1, 2, . . . , l in the matrix Σ can be directly calculated by the error covariance matrix of the state estimation in the ith local filter. However, the cross-covariance matrix is difficult to get. In our work, we use the result given by Gao et al. base on UT [17].
where χ Secondly, in MSIF, the optimal estimationX op theoretically lies in the closed interval: where norm2(•) denotes the 2-norm. If (50) is not maintained, the following degenerative methods (52) or (53) can be used:

The Proposed Method with Both Adaptivity and Robustness
In Section 2, we proposed the innovative algorithm to adapt Q when process-error exists. In this section, based on the hypothesis theory, the adopted judging criterion on process-error detection is given [4,17]. However, the drawbacks of the methods proposed by [4] and [17] are that the Q and R could not be estimated at the same time. Severe problems would be caused under the dilemma to decide which one should be adapted: Q, R, or both? To conquer this challenge, the RMNCE is employed to estimate the noise covariance of each sensor, then a decision is made whether the matrix R is suffering gross errors or not. Based on RMNCE, the decision is made easier by MSIF because the estimations of R are obtained with relatively high accuracy.
Briefly, in our proposed MSIF architecture, the matrix R of each sensor is estimated by RMNCE, and is denoted asR i k ; meanwhile, the process-error can be corrected by adapting Q if necessary. For the sake of simplicity, only one sensor of (1) is taken into consideration, so the superscript i inR i k could be dropped.

Robust R Estimaiton Based on RMNCE
As described in Section 2.2, the main advantage of RMNCE is that the estimate of variance is based only on measurements and hence can be immune to the state estimation error [27]. So,R k calculated by (17) can be considered as the benchmark to test whether the process error exists or not. Let R in be calculated by the following equation [28] The difference between R in andR k could be large if process error exists. The quotient is used as an indicator to detect the process error as where tr R k and tr(R in ) denote the trace ofR k and R in , respectively. When the quotient is around 1.0, it could be deduced with great probability that there is no process error. In this paper, a closed interval from 0.90 to 1.10 is used as the normal range. Otherwise, if the indicator is larger than 1.10 or smaller than 0.90, there is process error with great probability.

Remark 5.
If the diagonal elements differ by more than one order of magnitude, it is better to calculate the indicator separately. For example, the matricesR k and R in are aŝ Although the indicator calculated by (54) belongs to the closed interval, tr(R k ) tr(R in ) = 10.01 11 = 0.918 ∈ [0.90, 1.10], it is obvious that the second element in the diagonal differs 10 times in R in thanR k . Remark 6. The sliding window width for calculating R in andR k should be the same. In this paper, it is set to 50.
In our MSIF architecture, we also proposed the following matrix weights A R i as an alternative to (48). Consider a MSIF system with three sensors, the noise covariance matricesR i , i = 1, 2, 3 could be obtained by (13). Because the noise covariance matrices are usually diagonal, the following equation is employed to calculate the A R i where m is the dimension of Z i .

Remark 7.
In our proposed MSIF structure, all the sensors have the same measurement model, so all the A R i have the same dimension. If a MSIF system has different sensors, the variation of (56) should be derived specifically. However, the core idea is that the smallerR i (j, j) is, the greater weight A R i (j, j) is. A degenerated case is that if one or more elements of X are observed only by one sensor, assuming sensori, then where r denotes the number of elements in X that are only observed by sensor i. Hence, this is not a MSIF anymore, and this paper would not discuss this case any further.

Remark 8.
When there are l ≥ 2 sensors in a MSIF system, there are basically two options to utilize (13):

Option 1.
To calculate two matricesR k one time, so (13) and the relative equations will be run l − 2 times at least.

Option 2.
To calculate all matricesR k one time, then the classical least square method is used to solve the overdetermined equations, which are derived by the variation of (13). In this paper, Option 1 is adopted.

The Adaptive and Robust UKF Algorithm for MSIF (ARUKF-MSIF)
In this subsection, the complete scheme is given. The proposed Algorithm 1 ARUKF-MSIF aimed at target tracking can be implemented as follows. Step 1: State prediction through (5)-(7) for each sensor.
Step 3: Estimate matrix R i for each sensor through (12) to (17).
5.1 If (20) and (54)  The framework of the proposed ARUKF-MSIF methodology is shown in Figure 1. It has a two-layer fusion structure, the local layer and the globally optimal layer. Sensors 2021, 21, 5808 16 of 25

Simulations and Discussion
In this section, a set of numerical simulations of the radar tracking problems will be presented to illustrate the effectiveness of the proposed ARUKF-MSIF.

Process Model and Measurement Model
Similar to [1], consider the radar tracking system with three sensors. The differences between our model and [1] are threefold. Firstly, our model is more complicated and challenging: a two-dimension trajectory with time-varying process-error and/or measurement noise. In order to illustrate the effectiveness of our proposed method, the uncertainty of the system contains not only the measurement outliers/failures argued in [1] but also the process-error as in [15], which is usually the case in practice.
Secondly, all the sensors have the same measurement model in our simulation. With the development of hardware technology, more and more systems can afford redundant equipment to apply the multi-sensor information fusion into practical projects because the price of hardware is getting cheaper and cheaper. As proved by [6], if the sensors in decentralized structure have identical measurement matrices as in the centralized one, the two fusion methods are functionally equivalent. By taking advantage of the much smaller computational burden of the decentralized structure and the same accuracy as the centralized one, the simulation is arranged with three identical sensors whose noise characteristics are different.
Thirdly, the 2-D motion is taken into consideration as with [15,31]. The function of the radar can be demonstrated well enough; meanwhile, the model is not so complicated as the 3-D one. The simulated target trajectory is depicted in (Figure 2), and the true acceleration is shown in (Figure 3).
Consider that a target trajectory is in x-y plane. The position, velocity, and acceleration of the objective at time k are represented by the vector X k = x k , y k , .
x k , .. y k in the Cartesian coordinates. The dynamic equation for the target movement is as [15,31], where T denotes the sampling period. The initial state is X 0 = 1000 m, 5000 m, 10 m/s, 50 m/s, 2 m/s 2 , −4 m/s 2 , the process noise covariance matrix is The multi-sensor observation systems are composed by three radars with the same measurement model described as where r i k and ϕ i k denote the slant range and azimuth angle in polar coordinates measured by the ith radar, respectively. The variance matrices of V i k , i = 1, 2, 3 are different from each other. computational burden of the decentralized structure and the same accuracy as the centralized one, the simulation is arranged with three identical sensors whose noise characteristics are different.
Thirdly, the 2-D motion is taken into consideration as with [15,31]. The function of the radar can be demonstrated well enough; meanwhile, the model is not so complicated as the 3-D one. The simulated target trajectory is depicted in (Figure 2), and the true acceleration is shown in (Figure 3).  tralized one, the simulation is arranged with three identical sensors whose noise characteristics are different.
Thirdly, the 2-D motion is taken into consideration as with [15,31]. The function of the radar can be demonstrated well enough; meanwhile, the model is not so complicated as the 3-D one. The simulated target trajectory is depicted in (Figure 2), and the true acceleration is shown in (Figure 3).

Simulations
To demonstrate the effectiveness of our proposed method, the following three cases are designed specifically.

Case 1:
The noise covariance matrices R i = diag i 2 × 100, i 2 × 0.01 , i = 1, 2, 3 are known and the process noise covariance matrix Q varies over time. The process noise covariance matrix is assigned to be Q = diag[0.015, 0.015] during the epochs [200, 400].

Case 2:
The measurement noise covariance matrices R i are uncertain, and the process noise covariance matrix Q = diag[0.001, 0.001] is known. The measurement noise covariance matrix R 1 is assigned to be R 1 = 10 × diag 1 2 × 100, 1 2 × 0.01 during the epochs [200,400], and it is set to the same value as in Case 1 for the remaining periods.

Case 3:
Not only the measurement noise covariance matrices R 1 but also the process noise covariance matrix Q are uncertain. The changes in Case 1 and Case 2 are implemented simultaneously in this case.
Different from [15], in which the dynamic maneuver is simulated during a period when both the R i and the Q stay unchanged, in our simulations, all cases are designed to set the interval of the maneuver to the same one when R 1 and/or Q varies. So, the situations in our simulations are much harsher than the ones in [15]. The interacting multiple model (IMM) is employed by Ge et al. [15] to compare the algorithms; however, in this paper, we do not take the IMM method into consideration because its computational load is higher than the algorithm proposed by [15], and the complexity of our algorithm to adapt the Q is similar to [15]. Moreover, the Markov transition matrix used in IMM is obtained on the basis of statistics on the system evolution [33], so in order to employ the IMM, the transition matrix should be available. However, it could not be obtained precisely in practice.
For Case 1, the reason why only the matrix Q is uncertain is that we want to illustrate the effectiveness of our proposed method contrast to several algorithms which focus on the Q adaption. Except for the standard UKF, the AUKF algorithm proposed by [19], the adaptive fading UKF (AFUKF) method developed by [18], and the N-UKF given by [32] are taken into consideration together with ours to run the simulation.
For Case 2, to test the validity of our proposed MSIF structure with RMNCE, only R 1 changes during the specified intervals. The algorithms contained in the simulation are the improved Sage-Husa adaptive method in [26], the N-UKF, the standard UKF, and ours.
It is natural for Case 3 that the algorithms with both adaptability and robustness claimed by their inventors are taken into consideration except for the standard UKF and the improved Sage-Husa adaptive method. So, the N-UKF, the robust adaptive UKF proposed by [34], and our proposed method are applied to track the target.
For each case, R 1 of the sensor 1 is used by every single method. From 50 times of the Monte Carlo simulation in all cases, the root mean square error (RMSE) of each moment is obtained as where is the total times of simulation, x i k ,ŷ i k represents the filtering position of the target at time instant k in the ith simulation.
For the first case, the position tracking errors of the standard UKF, adaptive fading UKF, N-UKF, and our proposed method are shown in (Figure 4). To illustrate the performance of these algorithms clearly, we list the means and variance of the position tracking errors during the epochs [200, 600] and [601, 1000].
The epochs of Case 1 can be divided into three intervals: [0,199], [200,600], and [601, 1000]. The positioning errors during the first interval are close to each other because there is not process-modeling error. In the second interval, the maneuver of the target deteriorates the performance of the standard UKF. For the robust adaptive UKF, a modulation is made to prevent the algorithm divergence, i.e., a limitation is used for the adaptive factor. It can be seen clearly that all the methods except for the standard UKF can maintain their accuracy and robustness due to the Q adaption strategy used. From Table 1, the conclusion can be made that our proposed algorithm can resist the influence from both the time-varying process and the maneuvering motion models. The positioning errors estimated by our method remained the lowest among all the algorithms.
For Case 2, the improved Sage-Husa adaptive algorithm proposed by Zhang et al. is introduced to verify the performance in R estimation [26]. The position errors estimated by the algorithms mentioned above are shown in (Figure 5), the means and variances of the position tracking errors are listed in Table 2. The measurement noise variances estimated in the algorithms are shown in Figure 6.  introduced to verify the performance in R estimation [26]. The position errors estimated by the algorithms mentioned above are shown in (Figure 5), the means and variances of the position tracking errors are listed in Table 2. The measurement noise variances estimated in the algorithms are shown in Figure 6.    It can be seen from Figure 5 and Table 2 that the error of the improved Sage-Husa is much greater than other algorithms. This is because the innovation sequences used by this algorithm are contaminated by the maneuvering motion. The improved Sage-Husa method cannot maintain its accuracy or even diverges when process-error becomes nonnegligible, although it can overcome the time-varying noise covariance of the measurement. Similar to N-UKF, the robust adaptive UKF can recognize the mismatch between   It can be seen from Figure 5 and Table 2 that the error of the improved Sage-Husa is much greater than other algorithms. This is because the innovation sequences used by this algorithm are contaminated by the maneuvering motion. The improved Sage-Husa method cannot maintain its accuracy or even diverges when process-error becomes nonnegligible, although it can overcome the time-varying noise covariance of the measurement. Similar to N-UKF, the robust adaptive UKF can recognize the mismatch between the theoretical values and the calculated ones of Q and/or R , but it cannot provide the same accuracy as our proposed method because it cannot estimate the R independently. It can be seen from Figure 5 and Table 2 that the error of the improved Sage-Husa is much greater than other algorithms. This is because the innovation sequences used by this algorithm are contaminated by the maneuvering motion. The improved Sage-Husa method cannot maintain its accuracy or even diverges when process-error becomes non-negligible, although it can overcome the time-varying noise covariance of the measurement. Similar to N-UKF, the robust adaptive UKF can recognize the mismatch between the theoretical values and the calculated ones of Q and/or R, but it cannot provide the same accuracy as our proposed method because it cannot estimate the R independently.
From Figure 6, we can see that the improved Sage-Husa algorithm could estimate the measurement noise variance. However, this method requires a much longer time to catch up with the reference values; moreover, this method would give a much greater overestimation than ours. The N-UKF and the adaptive fading algorithm can estimate the measurement noise variance only when the variance becomes greater, and they fail when the variance becomes smaller. After the epoch 600, the estimated values cannot be restored to the reference values. The values estimated by the standard UKF are unchanged, because the standard UKF does not have the ability to estimate the measurement noise variance.
For Case 3, from Figure 7 and Table 3, we can see that the standard UKF takes about 200 epochs (from epochs [601, 800]) longer to obtain the accuracy than the N-UKF, the adaptive fading UKF, and our proposed method. The reason is that in order to maintain its accuracy, the standard UKF could only rely on the accurate measurement under the condition in Case 3, this correction procedure has a relative long lag. As to the improved Sage-Husa method, it could be influenced relatively significantly by the process-error more than other algorithms. The simulation result is similar as with Case 2, that the positioning errors estimated by the improved Sage-Husa algorithm are the greatest among the candidates. The N-UKF and the adaptive fading UKF have a similar accuracy because they can both resist the process-error and the time-varying measurement noise to some extent, but they cannot estimate the measurement noise variance as accurate as ours, as mentioned above in Case 2. to the reference values. The values estimated by the standard UKF are unchanged, because the standard UKF does not have the ability to estimate the measurement noise variance.
For Case 3, from Figure 7 and Table 3, we can see that the standard UKF takes about 200 epochs (from epochs [601, 800]) longer to obtain the accuracy than the N-UKF, the adaptive fading UKF, and our proposed method. The reason is that in order to maintain its accuracy, the standard UKF could only rely on the accurate measurement under the condition in Case 3, this correction procedure has a relative long lag. As to the improved Sage-Husa method, it could be influenced relatively significantly by the process-error more than other algorithms. The simulation result is similar as with Case 2, that the positioning errors estimated by the improved Sage-Husa algorithm are the greatest among the candidates. The N-UKF and the adaptive fading UKF have a similar accuracy because they can both resist the process-error and the time-varying measurement noise to some extent, but they cannot estimate the measurement noise variance as accurate as ours, as mentioned above in Case 2.

Discussion
Base on the decentralized MSIF architecture, we proposed the RAUKF-MSIF to tackle the problems for the nonlinear target tracking systems with time-varying noise covariance. From the results simulated in Case 1, we can see that if there are time-varying noise covariances, the algorithms can detect the process-error and adapt the Q except for the standard UKF.
In practice, the measurement noise covariance is also time-varying. So, in Case 2, we simulate the condition under which each sensors' measurement noise covariance is set to be 10 times greater during the specified interval than the others. From Figure 6, we can see that the estimation of R by our proposed scheme is more accurate than the others due to RMNCE and the fusion strategy described in Section 2.2 and 3.1. The performance of the improved Sage-Husa is not as good as ours because it could not resist the process-error.
The harshest environment of the simulations is presented in Case 3, in which, together with maneuvering, the matrix Q and R are changed simultaneously in the same interval. As analyzed in Section 4.2 for Case 3, our proposed method can obtain a relatively high ac-curacy due to both the process-error and the measurement noise variance can be estimated or adapted by our ARUKF-MSIF architecture.
Future work can be divided into two aspects. Firstly, we will focus on improving the numerical stability when adapting the matrix Q. Although the algorithms including ours can adapt the matrix Q, the adapted Q would be a negative definite or the adapted scalar factor would be so much greater than desired that some certain limited values must be set to prevent the algorithms from diverging.
Secondly, in addition to adapting the matrix Q, we would introduce the idea of IMM to estimate the dynamics of the target, while reducing the computational load of the IMM. Due to the existence of redundant measurement, we could improve the accuracy of the probability of the target moving from one model to another. Thus, a more precise model together with an adapted Q would improve the estimation accuracy further.

Conclusions
This paper addresses the target tracking problem with both time-varying process-error and time-varying measurement noise covariance by using a multi-sensor information filter structure based on an adaptive UKF and the RMNCE. The proposed ARUKF-MSIF employs the RMNCE to estimate the measurement noise variance of each sensor. Subsequently, it uses the theory of hypothesis testing to identify process model uncertainty; next, an innovative algorithm uses both the innovation sequences and the residual sequences to correct the process-error. Then, the algorithm calculates the weight matrices, which are used to fusion the results obtained by each individual sensor. The simulation results demonstrate that the proposed architecture can be robust against the process model uncertainty and error. Moreover, in the harshest circumstances as simulated in Case 3, the proposed method would still maintain the lowest position error among the algorithms. Further, although this paper uses the UKF as the filter, the core idea of the ARUKF-MSIF can be to broaden the EKF, CKF, and other variants of Kalman filter.