A Comprehensive Review of Multi-Source Data Fusion Processing Methods

1. Introduction

Based on Global Navigation Satellite Systems (GNSS), Inertial Navigation Systems (INS), Ultra-Wideband (UWB) technology, Bluetooth, Wireless Local Area Networks (WLAN), visual sensors, pseudolites (PL), and various other sensors, it is difficult to meet navigation performance requirements when used individually. In particular, GNSS, as a non-autonomous navigation system, has limitations in specific complex environments (such as urban canyons, tunnels, etc.), because its signals are prone to being blocked, interfered with, or shielded. An effective way to enhance the overall performance of a navigation system is to employ integrated navigation technology. This involves using two or more different types of navigation systems to measure and calculate the same navigation information, forming quantitative measurements. From these measurements, the errors of each navigation system are calculated and corrected. By employing a variety of technical means and methods, high accuracy and reliability of navigation and positioning services can be achieved in various scenarios, including seamless indoor and outdoor positioning, environments with electromagnetic interference, and underwater and underground environments. Therefore, multi-source data fusion positioning is based on the collaboration of multiple technology sources and adopts certain optimization criteria to achieve the best fusion positioning. The fusion method is the prerequisite and foundation for all-source navigation, and also serves as the key and core of integrated navigation systems.

Data fusion was initially introduced in an article titled Extension of the Probabilistic Data Association Filter in Multi-target Tracking published by the American system scientist Bar-Shalom. The article proposed the probabilistic data interconnection filter as the symbol of the formation of multi-source information fusion technology. The multi-source data fusion methods mainly adopt the switching method, the average weighted fusion method, and the adaptive weighted fusion method.

Based on the performance of different positioning sources, the optimal single positioning source is selected as the positioning means (Zou et al., 2013). However, ignoring other positioning sources is a waste and not the best choice. The average weighted fusion method does not take into account the different performances of different positioning sources but assigns the same weight to all positioning sources for fusion localization (Bhujle et al., 2016), which cannot achieve the optimal fusion effect. The adaptive weighted fusion method assigns different weights according to the characteristics of different fusion sources to achieve the best fusion positioning (Guo et al., 2018). The algorithms corresponding to these three methods mainly include optimal estimation algorithms, weighted fusion algorithms or adaptive weighted fusion algorithms (Yang et al., 2001; Yang et al., 2006; Zhang et al., 2016), Bayesian filters (BF), and variable-decibel Bayesian adaptive estimation (Sun et al., 2016; Li et al., 2023), Particle Filters (PF), Statistical decision theory (Berger et al., 2002), evidence theory (Zhang et al., 2017), fuzzy logic (Arikumar et al., 2016), etc. However, these algorithms all have specific preconditions and application scenarios, and it is necessary to establish an mathematical model between the observation information of the navigation source and the system state parameters.

In the field of dynamic positioning such as autonomous driving and vehicle navigation, the Kalman Filter (KF) has been widely used due to the introduction of physical motion models. However, KF is primarily designed for linear systems. For nonlinear systems, such as inertial navigation, the Extended KF (EKF) is suitable for weakly nonlinear objects because higher-order terms above the second order are discarded in the linearization process. To address the issue that the batch processing of the EKF random model requires storing a large amount of data, a recursive method for the random model has been proposed (Zhang et al., 2021), and time-domain non-local filtering data fusion algorithms have also been included (Cheng et al., 2017). The Unscented KF (UKF) retains the accuracy achieved by the third-order term of the Taylor series, making it suitable for nonlinear object estimation, although it involves relatively high computational demands. When both the system state and measurement noise are nonlinear, the PF can be used for nonlinear systems and systems with uncertain error models. However, the PF requires a probability density that closely approximates the real density. Otherwise, the filtering effect may be poor or even divergent. To address this, the Unscented Kalman Particle Filter (UPF) algorithm has been developed (Qin, 2020). However, both the PF and UFP methods face the issue of rapidly increasing computational load as the number of particles grows.

With the increasing demand from users for more comprehensive and intelligent navigation and positioning performance, filtering methods such as Factor Graphs (FG) and neural networks have been introduced. For example, FG algorithms have been extensively applied in single GNSS positioning, GNSS/INS integrated positioning, ambiguity resolution, and robust estimation (Watson et al., 2018; Wen et al., 2021; Gao et al., 2022; Zhang et al., 2024). To enhance positioning accuracy in urban environments, the FG algorithms have been optimized and improved (Lyu et al., 2021; Bai et al., 2022; Liao et al., 2023; Zhang et al., 2023). These studies have demonstrated that under certain conditions, FG algorithms exhibit higher computational accuracy and robustness compared to EKF. In 1965, Magill proposed the Multi-Model Estimation (MME) method, which enhances the adaptability of system models to real systems and changes in external environments under complex conditions, thereby improving the accuracy and stability of filtering estimates. The design of the model set, the selection of filters, estimation fusion, and the reset of filter inputs are all very important aspects. To enhance the high fault-tolerance capability of integrated navigation systems, Carlson introduced the Federated Filtering (FF) theory in 1988. This theory has been applied in indoor navigation, robotic navigation, and vision-language tasks (Carlson, 1988; Alsanhi et al., 2023). Existing Artificial Intelligence (AI) algorithms mainly include fuzzy control adaptive algorithms and neural network adaptive algorithms (Bian et al., 2004; Gao et al., 2007). For example, to address the impact of random disturbances on systems in underwater environments, RBF neural network-assisted FF has been employed for information fusion (Li et al., 2011). By establishing a black-box model with sufficiently accurate samples through offline training, the positioning accuracy and adaptability of the algorithm have been improved. To tackle the issues of high cost and susceptibility to weather conditions in existing high-precision satellite navigation for agricultural machinery, Yu et al. (2021) proposed a multi-sensor fusion automatic navigation method for farmland based on D-S-CNN. However, these AI algorithms require extensive training data, comprehensive pre-training of the system, and significant computational resources, and often struggle to ensure real-time performance, typically being used for post-processing.

Currently, scholars from various countries have conducted extensive research on integrated navigation systems. For instance, researchers from Linköping University in Sweden proposed a combined navigation system that integrates GPS, INS, and visible light vision assistance (Conte et al., 2009). This system utilizes the vision system and INS for positioning when GPS fails. Locata Corporation in Australia has integrated the Locata system with GPS, INS, vision systems, and Simultaneous Localization and Mapping (SLAM), achieving high-precision applications of the Locata system in both indoor and outdoor environments (Montillet et al., 2014). A communication and navigation fusion system has been applied for seamless positioning across wide-area indoor and outdoor spaces (Deng et al., 2013). A multi-frequency ground-penetrating radar data fusion system is used for working antennas in different frequency ranges (Coster et al., 2018), while multi-sensor data fusion is employed for analyzing airspeed measurement fault detection in drones (Guo et al., 2018). Additionally, an indoor mobile robot based on dead reckoning data fusion and fast response code detection (Nazemzadeh et al., 2017), and an IoT-based multi-sensor fusion strategy for analyzing occupancy sensing systems in smart buildings have been developed (Nesa et al., 2017). Systems that integrate vision, inertial navigation, and asynchronous optical tracking with Inertial Measurement Units (IMU) have also been implemented. Furthermore, several research teams have successfully developed open-source integrated navigation systems for use by academic or industrial technical personnel (Chen et al., 2021; Chi et al., 2023).

Although the aforementioned literature has conducted extensive research and testing on multi-source data fusion methods, fusion systems, and their applications. The theories and models of these methods have their specific applicable scenarios and conditions. Therefore, this paper summarizes the fundamental principles and mathematical models of multi-source data fusion methods, analyzes the advantages and disadvantages of different fusion approaches, and provides theoretical support and reference for the design, development, and application of fusion systems.

2. Multi-Source Data Fusion Processing Method

Fusion methods are primarily categorized into optimal estimation methods, filtering methods, MME, FG, FF, and other fusion approaches. The following sections will focus on elaborating the basic principles of these methods, their corresponding algorithmic models, and applicable scenarios.

2.1. Optimal Estimation Method

We estimate the parameter through randomly distributed observation vectors, specifically by finding a mapping function to compute the estimated value. Here, the system state is represented by the vector x, and the measurements corresponding to different navigation sources are denoted by y. The nonlinear and linear observation equations can be formulated as:

$$y=Ax+\epsilon$$

(1)

where $\mathit{A}$ is the observation matrix. $\mathit{\epsilon }$ includes the error term caused by random observation noise and linearization. The state of the system can be estimated by observation, and the estimated result is a vector $\widehat{\mathit{x}}$.

By employing different estimation criteria to address the problem of estimating unknown parameters, various estimation methods can be derived. Based on Equation (1) and the optimization of the criterion function, methods such as Least Squares Estimation (LSE), Minimum Variance Estimation (MVE), Maximum Likelihood Estimation (MLE), and Maximum A Posteriori Estimation (MAPE) can be formulated.

2.1.1. General LSE and Weighted LSE (WLSE)

The LSE is a parameter estimation algorithm proposed by the German mathematician Carl Friedrich Gauss in 1795, initially developed to determine planetary orbits (Koch, 1977). Regardless of whether the variable $x$ in the linear model (2) possesses prior statistical information or the random distribution followed by $\mathit{y}$, the LSE criterion employs a quadratic minimum.

$$J\left(\widehat{x}\right)={\left(y-Ax\right)}^{\mathrm{{\rm T}}}W\left(y-Ax\right)=\min$$

(2)

where $\mathit{W}$ is the appropriate positive definite weighted matrix, the above formula is the general LSE criterion when $\mathit{W}=\mathit{I}$ Then the estimation error and its variance are,

$$\left\{\begin{array}{l}{\widehat{x}}_{LS}={\left(A^{\mathrm{{\rm T}}}WA\right)}^{-1}{A}^{\mathrm{{\rm T}}}Wy\\ {\Sigma }_{\delta {\widehat{x}}_{LS}}={\left(A^{\mathrm{{\rm T}}}WA\right)}^{-1}{A}^{\mathrm{{\rm T}}}W{\Sigma }_{\epsilon }WA{\left(A^{\mathrm{{\rm T}}}WA\right)}^{-1}\end{array}\right.$$

(3)

It is easy to prove that when $W={\Sigma }_{\epsilon }^{-1}$，LSE has a minimum variance (Cui et al., 2001; Koch 1988a).

The most significant feature of this method lies in its algorithmic simplicity and independence from statistical information related to estimators and measurements. Previous researchers have implemented this algorithm in multi-source data fusion processing (Zhong et al., 2003; Jiao et al., 2010). However, the LSE exclusively utilizes measurement information for current state estimation, a characteristic that paradoxically constrains its applicability. Furthermore, while LSE optimality criterion guarantees the minimization of the total mean squared error in measurement estimates, it fails to ensure optimal estimation error for the estimator itself, consequently leading to suboptimal estimation accuracy.

2.1.2. MLE

Let the conditional probability density of $\mathit{y}$ with respect to $\mathit{x}$ be $f\left(\mathit{y}/\mathit{x}\right)$ and the MLE criterion be,

$${\left.f\left(y/x\right)\right|}_{x={\widehat{x}}_{ML}}=\max$$

(4)

This means that the $\mathit{y}$ has its maximum value at ${\left.f\left(y/x\right)\right|}_{x={\widehat{x}}_{ML}}=\max$, which is equivalent to,

$$𝜕f\left(y/x\right)/𝜕x{|}_{x={\widehat{x}}_{ML}}=0$$

(5)

Obviously, the solution of equation (5) is related to the $f\left(y/x\right)$. Different conditional probability density functions lead to different valuation formulas, therefore a universal formula cannot be derived. When the $f\left(y/x\right)$follows a normal distribution, that is,

$$f\left(y/x\right)={\displaystyle \frac{1}{{\left(2\pi \right)}^{n/2}{\left|{\mathsf{\Sigma }}_{y/x}\right|}^{1/2}}}\exp \left\{-{\displaystyle \frac{{\left(y-E\left(y/x\right)\right)}^{\mathrm{{\rm T}}}{\mathsf{\Sigma }}_{y/x}^{-1}\left(y-E\left(y/x\right)\right)}{2}}\right\}$$

(6)

According to the normal distribution condition expectation and variance formula (Koch 1988a; Cui et al., 2001),

$$\left\{\begin{array}{l}E\left(y/x\right)=E\left(y\right)+{\Sigma }_{yx}{\Sigma }_x^{-1}\left(x-E\left(x\right)\right)\\ {\Sigma }_{y/x}={\Sigma }_y-{\Sigma }_{yx}{\Sigma }_x^{-1}{\Sigma }_{xy}\end{array}\right.$$

(7)

From the above derivation, it can be concluded that MLE does not require prior distribution information about x. However, when its prior distribution information is known, the $f\left(y/x\right)$ should be derived strictly based on Equation (5).

2.1.3. MAPE

The estimation ${\widehat{x}}_{MA}$ of MAPE makes the $f\left(x/y\right)$ reach the maximum value under the condition $\widehat{x}={\widehat{x}}_{MA}$. The equivalent MAPE criteria are as follows after similar derivation (Cui et al., 2001):

$${\left.f\left(x/y\right)\right|}_{\widehat{x}={\widehat{x}}_{MA}}=\max$$

(8)

According to the conditional probability formula, it is easy to prove that (8) is equivalent to,

$${\left.f\left(x/y\right)\right|}_{\widehat{x}={\widehat{x}}_{MA}}=f\left(x/y\right)/f_2\left(y\right){|}_{\widehat{x}={\widehat{x}}_{MA}}=\max \iff f\left(x,y\right){|}_{\widehat{x}={\widehat{x}}_{MA}}=\max$$

(9)

That is, the MAPE uses the joint probability density function of the parameters and the observed data as the maximization criterion.

When both y and $x$ are normally distributed, the conditional probability density function is,

$$f\left(x/y\right)={\displaystyle \frac{1}{{\left(2\pi \right)}^{t/2}\left|{\Sigma }_{x/y}\right|}}\exp \left\{-{\displaystyle \frac{{\left(x-E\left(x/y\right)\right)}^{\mathrm{{\rm T}}}{\Sigma }_{x/y}^{-1}\left(x-E\left(x/y\right)\right)}{2}}\right\}$$

(10)

Differentiate Equation (9) and set the derivative to zero to obtain the variance and estimation error of the MAPE:

$$\left\{\begin{array}{l}{\widehat{x}}_{MA}=E\left(x/y\right)\\ {\Sigma }_{\delta {\widehat{x}}_{MA}}={\Sigma }_{x/y}\end{array}\right.$$

(11)

When x is uncorrelated with y, the estimated result is equal to the prior information of the parameter.

Compared with MLE, MAPE places greater emphasis on the prior information of parameters. Therefore, MAPE becomes meaningless when parameters are non-random quantities. In practice, MAPE can be understood as a modification of parameters through observational data, where the extent of this modification is determined by the variance of the observations and their correlation with the parameters.

2.1.4. MVE

The MVE is an optimization criterion that minimizes the variance of the estimation error to estimate ${\widehat{x}}_{MV}$.

$$J\left(\widehat{x}\right)={\left.E\left({\left(x-{\widehat{x}}_{MV}\right)}^{\mathrm{{\rm T}}}\left(x-{\widehat{x}}_{MV}\right)\right)\right|}_{\widehat{x}={\widehat{x}}_{MV}}={\displaystyle {\int }_{-\infty }^{\infty }{\displaystyle {\int }_{-\infty }^{\infty }{\left(x-{\widehat{x}}_{MV}\right)}^{\mathrm{{\rm T}}}\left(x-{\widehat{x}}_{MV}\right)}}f\left(x,y\right)dxdy=\min$$

(12)

According to the conditional and edge probability density functions, equation (12) is equivalent to

$$h\left({\widehat{x}}_{MV}\right)={\displaystyle {\int }_{-\infty }^{\infty }{\left(x-{\widehat{x}}_{MV}\right)}^{\mathrm{{\rm T}}}\left(x-{\widehat{x}}_{MV}\right)}dx=\min$$

(13)

The MVE and its estimation error are:

$$\left\{\begin{array}{l}{\widehat{x}}_{MV}=E\left(x/y\right)\\ {\mathsf{\Sigma }}_{\delta {\widehat{x}}_{MV}}={\displaystyle {\int }_{-\infty }^{\infty }{\mathsf{\Sigma }}_{x/y}}{f}_2\left(y\right)dy\end{array}\right.$$

(14)

where $f_2\left(y\right)$ is the edge probability density function of the observation. When both the x and the y are normally distributed, the MVE is completely equivalent to the MAPE.

2.1.5. LMVE

MLE, MAPE, and MVE all require the joint probability density or conditional probability density function of $\mathit{y}$ and $\mathit{x}$. Their estimation formulas are related to distribution information and are not necessarily linear combinations of the observed values. However, LSE does not require the distribution information of the parameters and observed values, the parameter estimation is a linear combination of the observed values and linear estimation. As the name suggests, The LMVE is also a linear estimation method, which does not require the distribution information of the observed values and parameters, but only their statistics.

Let LMVE be estimated as

$${\widehat{x}}_{LMV}=\alpha +\beta y$$

(15)

where $\alpha$ is a constant vector and $\beta$ is a constant matrix, then the estimation error and its variance are:

$$\left\{\begin{array}{l}\delta {\widehat{x}}_{LMV}=x-{\widehat{x}}_{MV}=x-\alpha -\beta y\\ {\Sigma }_{\delta {\widehat{x}}_{LMV}}={\Sigma }_x+\beta {\Sigma }_y{\beta }^{\mathrm{{\rm T}}}-{\Sigma }_{xy}{\beta }^{\mathrm{{\rm T}}}-\beta {\Sigma }_{xy}\end{array}\right.$$

(16)

It is worth noting that LMVE takes the minimum mean square error ${\Sigma }_{\delta {\widehat{x}}_{LMV}}$ as the estimation criterion.

$$\begin{array}{cc}{M}_{\delta {\widehat{x}}_{LMV}}\hfill & =E\left(\delta {\widehat{x}}_{LMV}^{\mathrm{{\rm T}}}\delta {\widehat{x}}_{LMV}\right)=E\left(\delta {\widehat{x}}_{LMV}\right)E{\left(\delta {\widehat{x}}_{LMV}\right)}^{\mathrm{{\rm T}}}+{\Sigma }_{\delta {\widehat{x}}_{LMV}}\hfill \\ \hfill & =E\left(\delta {\widehat{x}}_{LMV}\right)E{\left(\delta {\widehat{x}}_{LMV}\right)}^{\mathrm{{\rm T}}}+\left(\beta -{\Sigma }_{xy}{\Sigma }_y^{-1}\right){\Sigma }_y\left(\beta -{\Sigma }_{xy}{\Sigma }_y^{-1}\right)+{\Sigma }_x-{\Sigma }_{xy}{\Sigma }_y^{-1}{\Sigma }_{yx}\hfill \end{array}$$

(17)

Then,

$$M_{\delta {\widehat{x}}_{LMV}}=\min \iff \left\{\begin{array}{l}E\left(\delta {\widehat{x}}_{LMV}\right)=0\\ \beta ={\Sigma }_{xy}{\Sigma }_y^{-1}\end{array}\right.\iff \left\{\begin{array}{l}E\left(\delta {\widehat{x}}_{LMV}\right)=0\\ \mathrm{var}\left(\delta {\widehat{x}}_{LMV}\right)=\min \end{array}\right.\iff \mathrm{var}\left(\delta {\widehat{x}}_{LMV}\right)=\min$$

(18)

Therefore, LMVE is unbiased estimation, and the variance of estimation error is minimal under the premise of unbiased estimation. However, the minimum variance of the estimation error is only a necessary condition for LMVE. By inserting equation (18) into equation (15) and (16), The estimation error and variance of LMVE are:

$$\left\{\begin{array}{c}{\widehat{x}}_{LMV}=E\left(x\right)+{\Sigma }_{xy}{\Sigma }_y^{-1}\left(y-E\left(y\right)\right)={\mu }_x+{\Sigma }_{xy}{\Sigma }_y^{-1}\left(y-E\left(y\right)\right)\\ {\Sigma }_{{\widehat{x}}_{LMV}}={\Sigma }_{xy}{\Sigma }_y^{-1}{\Sigma }_{yx}\end{array}\right.$$

(19)

It is noteworthy that equation (19) is derived under the premise of knowing the statistical properties (namely, the mathematical expectation and variance) of y and x, and is independent of their distribution information. Therefore, when both x and y follow a normal distribution, the LMVE is equivalent to the MVE and MAPE.

2.1.6. Comparison of Several Different Optimal Estimation Methods

When x and y are normally distributed, Table 1 presents a table of parameter estimates based on different criteria.

Table 2 presents the advantages and disadvantages of several optimal estimation methods.

2.2. Filtering Algorithm

The optimization method utilizes the observations from the current epoch to estimate the system state, hence the localization results are significantly influenced by the current observation errors. In 1960, R. E. Kalman first proposed the KF, which employs a recursive approach to avoid the accumulation of large amounts of data and designs the filter using the state-space method in the time domain, making KF suitable for the estimation of multi-dimensional random processes. Depending on the differences in system equations and state equations, Filtering algorithm is primarily divided into the following algorithms.

2.2.1. Standard KF

Assuming that both the state motion and observation models are linear, the estimated state $x_k$ at epoch $t_k$ is driven by the system noise sequence $w_{k-1}$, and the driving mechanism is described by the following equation of state:

$$x_k={\Phi }_{k,k-1}{x}_{k-1}+{\Gamma }_{k-1}{w}_{k-1}$$

(20)

where ${\Phi }_{k,k-1}$ is a one-step transfer matrix from time $t_{k-1}$ to time $t_k$. ${\Gamma }_{k-1}$ is the system noise-driving term. Let the observation of $x_k$ be $z_k$ , $z_k$ and $x_k$ satisfy the linear relationship, and the observation equation is:

$$z_k=H_k{x}_k+v_k$$

(21)

where $H_k$ is the observation matrix, $v_k$ is the observation noise, $w_k$ and $v_k$ meet at the same time:

$$E\left(w_k\right)=0,\mathrm{cov}\left(w_k,w_j\right)=E\left(w_k{w}_j^{\mathrm{{\rm T}}}\right)=Q_k{\delta }_{kj}$$

(22)

$$E\left(v_k\right)=0,\mathrm{cov}\left(v_k,v_j\right)=E\left(v_k{v}_j^{\mathrm{{\rm T}}}\right)=R_k{\delta }_{kj}$$

(23)

$$\mathrm{cov}\left(w_k,v_j\right)=E\left(w_k{v}_j^{\mathrm{{\rm T}}}\right)=0$$

(24)

where $Q_k$ is the variance matrix of the state noise sequence. $R_k$ is the variance matrix of the observation noise sequence.

If the $Q_k$ is not negative definite, the $R_k$ is positive definite, then the estimate ${\widehat{x}}_k$ is solved according to the KF equation. The calculation steps of KF include a filter gain calculation loop (time update process) and a filter calculation process (measurement update process) two calculation loops. As shown in Figure 1. Areas with green shading represent the estimates and their covariance matrix.

As can be seen from Figure 1, the KF must provide initial valuesand. Practice has proven that the recursive computation method is the greatest advantage of KF. This algorithm can be executed by a computer and does not require the storage of large amounts of measurement data over time, which is why the KF has been widely adopted in engineering. KF clearly states that the system driving noise and measurement noise must be white noise. However, in reality, these two types of noise in some systems are colored noise. Consequently, the KF under colored noise conditions has been proposed, mainly including the KF when the system noise is colored and the measurement noise is white, and the KF when the system noise is white and the measurement noise is colored (Jazwinski, 1970).

The estimated $P_{k|k+1}$ calculated according to the filtering equations tends to zero or to a steady-state value as the number of observation epochs k increases. When the deviation between the estimated value and the true value becomes increasingly large, the filter gradually loses its estimation capability. This phenomenon is known as filter divergence. To mitigate filter divergence, numerous scholars have proposed various methods, including information filtering (Dan, 2006), square root filtering (Bierman et al., 1987), UUDT decomposition filtering (Li et al., 2009), adaptive filtering (Mechra et al., 1972), suboptimal filtering (Athans et al., 1968), and H∞ filtering (Dan, 2006). Each of these methods has its own advantages, disadvantages, and applicable scopes, and further details can be found in the relevant literature.

2.2.2. Extended Kalman Filter (EKF)

The mathematical model of the system is linear in the KF, which means it is only applicable when both the system equations and observation equations are linear. However, in engineering applications, the mathematical models of systems are often nonlinear, such as in the INS of aircraft and ships, satellite navigation, and industrial control systems, making KF unsuitable for direct application. In such cases, the EKF can be used to linearize the nonlinear system. EKF utilizes Taylor series expansion to linearize the nonlinear system and transforms the original nonlinear system model into state and observation equations.

The nonlinear equation of state and observation equation corresponding to EKF are:

$$\left\{\begin{array}{l}\dot{x}\left(k\right)=f\left(x\left(k\right)\right)+g\left(w\left(k\right)\right)\\ z\left(k\right)=h\left(x\left(k\right)\right)+v\left(k\right)\end{array}\right.$$

(25)

where $f\left(x\left(k\right)\right)$ and $g\left(w\left(k\right)\right)$ are nonlinear functions of the state vector $x\left(k\right)$ and system noise vector $w\left(k\right)$ respectively. $h\left(x\left(k\right)\right)$ is the nonlinear function of the design matrix.

The two equations in formula (25) are linearized as follows:

$$\left\{\begin{array}{l}\delta x_k=F_{k-1}\delta x_{k-1}+G_{k-1}{w}_{k-1}\\ \delta z_k=H_k\delta x_k+v_k\end{array}\right.$$

(26)

After linearization, the solution can be solved by using KF.

The EKF is suitable for weakly nonlinear systems because its linearization process retains only the first-order terms while discarding higher-order terms above the second-order. To address the issue of the EKF stochastic model batch processing requiring the storage of large amounts of data, the algorithm proposes a recursive approach for the stochastic model (Zhang et al., 2021). Additionally, it includes a time-domain non-local filtering data fusion algorithm (Cheng et al., 2017).

2.2.3. Unscented Kalman Filter (UKF)

The EKF discards higher-order terms above the second order and retains only the linear terms, thus the EKF algorithm is only suitable for the estimation of weakly nonlinear objects. When the nonlinearity of the estimated object is stronger, the resulting estimation error will also be larger, and may even cause the filter to diverge. The UKF is an effective method for solving nonlinear system problems. In 1995, S.J. Julier and J.K. Uhlmann proposed the UKF algorithm to address the filtering problem of strongly nonlinear objects, which was subsequently further refined by E.A. Weiss and R.V. Merwe (Simon, 2001; Dan, 2014).

The core of UKF is to use Untrace Transformation (UT) to determine the mapping relationship between variables, which is equivalent to retaining the accuracy achieved by the third-order term of the Taylor series, so it is suitable for nonlinear object estimation.

Using discrete nonlinear systems, i.e.:

$$\left\{\begin{array}{l}{x}_k=f\left(x_{k-1}\right)+w_k\\ z_k=h\left(x_k\right)+v_k\end{array}\right.$$

(27)

where $x_k$ is the state vector. $w_k$ is the system noise vector. $z_k$ is the observed vector. $v_k$ is the is the observed noise vector.

A series of sampling points are selected near $x_{k-1}$, The mean and covariance of these sample points are ${\widehat{x}}_{k-1}$ and $P_{k-1}$, respectively. The corresponding transform sampling points are generated through the nonlinear system. The predicted mean and covariance can be obtained by calculating these transform sampling points.

Let the state variable be n-dimensional, then 2n+1 sampling points and their weights are:

$$\left\{\begin{array}{l}{\xi }_0={\widehat{x}}_{k-1}\\ {\xi }_i={\widehat{x}}_{k-1}+\sqrt{n+\tau }{\left(P_{k-1}\right)}_i\\ {\xi }_{i+n}={\widehat{x}}_{k-1}-\sqrt{n+\tau }{\left(P_{k-1}\right)}_i\end{array}\right.i=1,2,\cdots ,n$$

(28)

and the weight corresponding to ${\xi }_{i,k-1}\left(i=0,1,\cdots ,2n\right)$ is:

$$W_i^m=W_i^c=\left\{\begin{array}{c}\kappa /\left(n+\kappa \right)i=0\\ 1/2\left(n+\kappa \right)i\ne 0\end{array}\right.$$

(29)

where $\kappa$ is the proportional coefficient, which can be used to adjust the distance between sigma points and ${\widehat{x}}_{k-1}$, and this coefficient only affects the higher-order matrix deviations above the second order. ${\left(\sqrt{P_{k-1}}\right)}_i$ is the i-th row or column of the square root matrix. There are different ways to determine $W_i^m$, and some literature has made some improvements, such as the UKF algorithm under additive noise cases and non-additive noise conditions (Liu et al., 2010). The specific process of the UKF algorithm can be referred to in the relevant literature. Practice proves that UKF is suitable for nonlinear object estimation, but the computation is relatively large.

Due to the strong advantages of the UKF in handling nonlinear systems, it has also been applied in multi-source data fusion. Examples include autonomous multi-level positioning based on smartphone-integrated sensors and pedestrian indoor networks (Shi et al., 2023), as well as the use of UKF in mobile mapping applications based on low-cost GNSS/IMU under demanding conditions (Cahyadi et al., 2024).

2.2.4. Particle Filter (PF)

The PF was initially proposed by Metropolis and Wiener as early as 1940 (Wiener, 1956). PF is a method that approximates the probability density function by propagating a set of random samples in the state space. It replaces the integral operation with the sample mean, thereby obtaining the minimum variance estimate of the system state. These samples are referred to as particles in the image, hence the method is known as PF. Here, the samples are the particles, and when the number of samples N→∞, it can approximate any form of probability density distribution.

The PF directly calculates the conditional mean, which is the minimum variance estimate, based on the probability density. This probability density is approximated by the EKF or the UKF. The estimate is determined by the weighted average of sample values (particles) from multiple different distributions. Each particle computation values (particles) from multiple different distributions. Each particle computation requires a complete EKF or UKF calculation. Therefore, PF is suitable for estimation under nonlinear system and measurement conditions, offering higher estimation accuracy than using EKF or UKF alone, but with a significantly higher computational load compared to EKF and UKF.

PF has many different names, For example, Sequential Importance Sampling (SIS) (Doucet et al., 2001), Bootstrap Filtering (Gray et al., 1993), Condensation Algorithm (Isaard et al., 1996; MacCormick et al., 1999), Interacting Particle Approximations (Moral, 1998), Monte-Carlo Filtering (Kitagawa, 1996), Sequential Monte-Carlo filtering (Crisan et al., 2002; Andrieu et al., 2004).

According to relevant literature, the general procedure for executing PF is as follows:

Step 1: Initial value determination

The initial particle value ${\chi }_0^{\left(i\right)}\left(i=1,2,\cdots N\right)$ is generated according to the prior probability density $p\left(x_0\right)$ of the initial state.

For k=1, 2, 3, ... Execute,

Step 2: Select the recommended probability density $q\left[x_k/\left(x_0^k\left(i\right),z_0^k\right)\right]$. According to this recommended density, a particle ${\chi }_k^{\left(i\right)}$ at time k, i=1,2,... N is generated, as a secondary sample of the original particle.

Calculated weight coefficient,

$$w_k^{\left(i\right)}=w_{k-1}^{\left(i\right)}{\displaystyle \frac{p\left[z_k/{\chi }_k^{\left(i\right)}\right]p\left[{\chi }_k^{\left(i\right)}/{\chi }_{k-1}^{\left(i\right)}\right]}{q\left[{\chi }_k^{\left(i\right)}/\left({\chi }_0^{\left(i\right)},z_0^k\right)\right]}}$$

(30)

$$w_0^{\left(i\right)}=p\left[{\chi }_0^{\left(i\right)}\right]$$

(31)

$${\tilde{w}}_k^{\left(i\right)}=w_k^{\left(i\right)}\left[{\displaystyle \sum _{j=1}^N{w}_k^{\left(i\right)}}\right],i=1,2,\cdots n$$

(32)

Step 3: the PF algorithm, such as the SIR method or residual secondary sampling method, is used to perform secondary sampling on the original particle ${\chi }_k^{\left(i\right)}\left(i=1,2,\cdots n\right)$ and generate secondary sampling to update the particle ${\chi }_k^{+\left(i\right)}\left(i=1,2,\cdots n\right)$.

Step 4: Calculating the filtering value based on the secondary sampling particles,

$${\widehat{x}}_k={\displaystyle \frac{1}{n}}{\displaystyle \sum _{i=1}^n{\chi }_k^{+\left(j\right)}}$$

(33)

The PF is mainly used for data fusion and integrity monitoring (Wang et al., 2014 & 2016; Yun et al., 2017; Bi et al., 2019). Therefore, when both system and measurement noises are nonlinear, PF can be applied to nonlinear systems as well as systems with uncertain error models.

2.2.5. UKF-Based Particle Filter (UPF)

The core of PF is to select the proposal probability density. The closer the proposal density is to the true density, the better the filtering effect will be. Conversely, if there is a significant difference between the proposal density and the true density, the filtering effect will deteriorate, and divergence may even occur. If PF is combined with UKF, and the proposal density is determined by UKF, the problem of particle degeneracy can be resolved. When updating the particles, the latest posterior information can be obtained, which is beneficial for regions with a high particle likelihood ratio. The method that combines PF with UKF is called UPF. However, both PF and UPF have a common problem, which is that as the number of particles increases, the computational complexity will sharply increase.

2.2.6. Federated Filtering (FF)

The filtering methods introduced above belong to Centralized KF (CKF). CKF has problems such as high state dimension, heavy computational burden, and poor fault-tolerance. Another filtering method is Decentralized KF (DKF). DKF has been developed for over 20 years. As early as 1971, Pearson (1971) proposed the concept of dynamic decomposition and a two-level structure for state estimation. Subsequently, Speyer (1979), Willsky (1982), Kerr (1987), and Carlson (1988) made contributions to DKF techniques. Among the many DKF methods, the FF proposed by Carlson has been valued for its design flexibility, low computational load, and good fault-tolerance. Now, the FF has been selected as the basic algorithm for the U.S. Air Force's fault-tolerant navigation system, the Common Kalman Filter program (Loomis et al., 1988).

The FF proposed by Carlson is designed to address the following issues:

(1) The filter should have good fault-tolerance. When one or several navigation subsystems fail, it should be able to easily detect and isolate the faults and quickly recombine the remaining normal navigation subsystems (reconfiguration) to continue providing the required filtering solution.

(2) The filtering accuracy should be high.

(3) The fusion algorithm from local filtering to global filtering should be simple, with low computational load and minimal data communication, to facilitate real-time implementation of the algorithm.

The FF is a two-level filter, as shown in Figure 2. The Carlson FF introduces a master filter. Since the master filter does not accept measurement inputs, it only has time updates and no measurement updates. Additionally, a feedback control switch from the master filter to the sub-filters is added. The entire filtering system consists of $\overline{N}=\left(N+1\right)$ filters, with N sub-filters providing local estimates (${\widehat{X}}_k^{\left(ci\right)}$ and $P_k^{\left(ci\right)}$) to the master filter. These local estimates are optimally combined with the master filter estimates (${\widehat{X}}_k^{\left(m\right)}$ and $P_k^{\left(m\right)}$) to obtain the global estimates (${\widehat{X}}_k^g$ and $P_k^g$).

If there are N local state estimates ${\widehat{X}}_1$, ${\widehat{X}}_2$,…, ${\widehat{X}}_N$ and their corresponding covariance matrices $P_{11}$, $P_{22}$, …,​ $P_{NN}$, and the local estimates are uncorrelated with each other, i.e., $P_{ij}=0\left(i\ne j\right)$, then the global optimal estimate can be expressed as:

$${\widehat{X}}_g=P_g{\displaystyle \sum _{i=1}^N{P}_{ii}^{-1}}{\widehat{X}}_i$$

(34)

$$P_g={\left({\displaystyle \sum _{i=1}^N{P}_{ii}^{-1}}\right)}^{-1}$$

(35)

Let the FF have N (N > 2) sub-filters, and the output of sub-filter i is,

$${\widehat{X}}_i=X+{\tilde{X}}_{i}i=1,2,\mathrm{...}N$$

(36)

where X represents the common state of all sub-filters, with a dimension of n. ${\tilde{X}}_i$​ denotes the estimation error of the $i$-th sub-filter. If the sub-filter operates normally, ${\tilde{X}}_i$ is white noise.

The measurement equation for X can be formulated based on the outputs of the N sub-filters.

$$Z=HX+V$$

(37)

Thus,

$$Z=\left[\begin{array}{c}{\widehat{X}}_1\\ {\widehat{X}}_2\\ ⋮\\ {\widehat{X}}_N\end{array}\right],H=\left[\begin{array}{c}{I}_{n\times n}\\ I_{n\times n}\\ ⋮\\ I_{n\times n}\end{array}\right],V=\left[\begin{array}{c}{\tilde{X}}_1\\ {\tilde{X}}_2\\ ⋮\\ {\tilde{X}}_N\end{array}\right]$$

(38)

Assuming that each sub-filter operates normally and the estimation errors are mutually uncorrelated, we have，

$$\left\{\begin{array}{l}E\left(V\right)=0\\ R=E\left(V{V}^{\mathrm{{\rm T}}}\right)=diag\left(\begin{array}{cccc}{P}_{11}& P_{22}& \cdots & P_{NN}\end{array}\right)\end{array}\right.$$

(39)

Among them, $P_{ii}=E\left({\tilde{X}}_i{\tilde{X}}_i^{\mathrm{{\rm T}}}\right)$ is the covariance matrix of the estimation error of sub- filter $i$.

According to the literature, the Markov estimate of the common state X is,

$$\left\{\begin{array}{l}{P}_g={\left(H^{\mathrm{{\rm T}}}{R}^{-1}H\right)}^{-1}={\left({\displaystyle \sum _{i=1}^N{P}_{ii}^{-1}}\right)}^{-1}\\ {\widehat{X}}_g={\left(H^{\mathrm{{\rm T}}}{R}^{-1}H\right)}^{-1}{H}^{\mathrm{{\rm T}}}{R}^{-1}Z={\left({\displaystyle \sum _{i=1}^N{P}_{ii}^{-1}}\right)}^{-1}{\displaystyle \sum _{i=1}^N{P}_{ii}^{-1}}{\widehat{X}}_i=P_g{\displaystyle \sum _{i=1}^N{P}_{ii}^{-1}}{\widehat{X}}_i\end{array}\right.$$

(40)

The physical meaning of the above result is quite evident. If the estimation accuracy of ${\tilde{X}}_i$ is poor, that is, $P_{ii}$ is large, then its contribution to the global estimate is relatively small. The above discussion pertains to the fusion algorithm when the estimates of each sub-filter are uncorrelated. For the fusion algorithm when the estimates of each sub-filter are correlated, one can refer to relevant literature.

2.2.7. Comparison of Different Filtering Method

To more clearly illustrate the performance of the various filtering methods introduced above, Table 3 provides a summary of the information.

2.3. Multiple Model Estimation (MME)

MME can estimate the system state through weighted first-order filter estimates with different parameter values, thereby achieving the goal of adapting to unknown or uncertain system parameters (Magill et al., 1965). In 1970, Ackerson first applied MME to jump environments, arguing that the system pattern is a finite-state Markov chain that can switch between different patterns. Since then, MME methods have been widely used in many fields under various names (Li, 1996), such as multi-model adaptive estimation, parallel processing algorithm, filter bank method, segmented filter, and improved Gaussian sum filter.

Any real system has different degrees of uncertainty, which is sometimes manifested inside the system and sometimes manifested outside the system. From the inside of the system, the mathematical model structure and parameters describing the controlled object cannot be fully known in advance by the designer. As for the influence of the external environment on the system, it can be expressed as many disturbances. These disturbances are often unpredictable and can be either deterministic or random. In addition, some measurement noise enters the system from different measurement feedbacks in the same way as disturbances, and the statistical properties of these random disturbance noises are usually unknown. Therefore, under the condition that the mathematical models of the controlled object and the disturbances are not fully determined, the control sequence is designed to make the specified performance index as close to and as optimal as possible. Essentially, an adaptive control system is an intrinsically nonlinear system, which is very difficult to analyze. MME theory employs N linear stochastic control systems to solve the nonlinear problem of adaptive control, which will improve the adaptability of the system models to real systems and external environmental changes, and enhance the accuracy and reliability of filter estimation.

Assume the linear stochastic system model is:

$$\left\{\begin{array}{l}x\left(k+1\right)=\Phi \left(\theta \right)x\left(k\right)+B\left(\theta \right)u\left(k\right)+\Gamma \left(\theta \right)w\left(k\right)\\ z\left(k\right)=H\left(\theta \right)x\left(k\right)+v\left(k\right)\end{array}\right.$$

(41)

where $x\left(k\right)$ is the state vector of the system, $z\left(k\right)$ is the output vector of the system, and $\Phi \left(\theta \right)$, $B\left(\theta \right)$ , $\Gamma \left(\theta \right)$, $H\left(\theta \right)$ are the system matrices of appropriate dimensions. $u\left(k\right)$ and $w\left(k\right)$ are sequences of zero-mean white noise vectors with dimensions n and m, respectively, and their covariance matrices are $Q\left(\theta \right)$ and $R\left(\theta \right)$.

The nonlinear system is linearized at N operating points, allowing the original nonlinear system to be approximated by N sets of linear equations. The parameter $\theta \in \left\{{\theta }_1,{\theta }_2,\cdots {\theta }_N\right\}$ can take N discrete values, thereby forming a combination of N linear systems that serve as an approximation of the nonlinear system model.

By representing any value of the parameter $\theta$ as ${\theta }_i$, the system matrix in equation (41) is redefined as follows:

$$\left\{\begin{array}{l}\Phi \left({\theta }_i\right)={\Phi }_i,B\left({\theta }_i\right)=B_i\\ \Gamma \left({\theta }_i\right)={\Gamma }_i，H\left({\theta }_i\right)=H_i\left(i=1,2,\cdots N\right)\end{array}\right.$$

(42)

Based on the aforementioned notation, N discrete random linear systems can be characterized as,

$$\begin{array}{l}x{}_i\left(k+1\right)=\Phi {}_{i}x\left(k\right)+B{}_{i}u\left(k\right)+\Gamma {}_{i}w\left(k\right)\\ z{}_i\left(k\right)=H{}_{i}x\left(k\right)+v{}_i\left(k\right)\end{array}$$

(43)

Figure 3 illustrates the fundamental principle of the MME method, wherein N discrete values constitute N distinct systems.

As illustrated in Figure 3, a parallel array of filters is implemented to accommodate different operational modes of the stochastic hybrid system. Each filter processes both the control input and measurement data from the system, generating output residuals and state estimates based on individual models. The system incorporates model probability design for each corresponding filter, and the comprehensive state estimation is derived from the weighted average of all filter state estimates. Table 4 summarizes the performance characteristics of the MME method.

To date, the MME method has been extensively utilized in diverse fields, including integrated navigation data fusion, self-calibration, and target tracking (Luan et al., 2021; Tian et al., 2022; Yang et al., 2024).

2.4. Factor Graph (FG) Methods

In the filtering method mentioned earlier, the state of the system at the current moment is only related to the observations at the current moment and the navigation state at the previous moment. However, in practical applications, some observations may be delayed, and certain position solutions need to be realized over a period before and after the observations. This cannot be described solely by the state information of the current and previous moments. For example, the EKF converts historical observations into prior information of the current state through the state equation for propagation, and the state linearization points corresponding to the historical observations after propagation are fixed. When there is an undetected blunder in the historical observations, the linearization point error is large, which can easily lead to the prior information being contaminated, thereby affecting positioning accuracy. In contrast, the FG can fully utilize historical observations by iteratively updating the linearization points, mining the constraint information of the observations in the temporal dimension, and suppressing the influence of blunders.

FG decomposes a complex global function involving multiple variables into the product of several simpler local functions, thereby constructing a bidirectional graph structure. When dealing with global functions that involve numerous variables, the conventional approach is to break down the given function into its constituent factors, which serve as local functions. These local functions are then combined multiplicatively to represent the original global function. The FG is a bipartite graph model that elegantly captures this factorization process. It typically consists of variable nodes, factor nodes, and edges connecting them. By leveraging this structure, FG can effectively decompose multivariable functions into products of local functions, facilitating more efficient computation and analysis.

$$g\left(x_1,x_2,\cdots x_n\right)={\displaystyle \prod _{j\in J}{f}_j\left(X_j\right)}$$

(44)

where $x_n$ denotes a variable node, the local function $f_j$ denotes a factor node, denotes the independent variable associated with the local function $f_j\left(\cdot \right)$, and subsequently, the factor node is connected to the corresponding variable node via an edge that represents their mutual relationship.

Let $g\left(x_1,x_2,x_3,x_4,x_5\right)$ be a global function involving five variables. It can be decomposed into the product of four factors, which can be expressed as follows:

$$g\left(x_1,x_2,x_3,x_4,x_5\right)=f_1\left(x_1,x_3\right)f_2\left(x_1,x_2,x_4\right)f_3\left(x_3,x_4\right)f_4\left(x_4,x_5\right)$$

(45)

The corresponding FG structure is shown in Figure 4.

Figure 5 establishes a factor graph-based multi-sensor fusion framework, incorporating sensors including IMU, GPS, barometric altimeter, optical flow sensor, magnetic heading sensor, and star tracker.

In Figure 5, the large circles represent variable nodes, and the black solid small squares represent factor nodes. The variables of the functions associated with the factor nodes include the variable nodes connected to them. Each factor has a corresponding error function. By adjusting the x to minimize the error of the FG, the optimal estimate can be obtained. The formula is as follows:

$$\widehat{x}=\underset{\mathrm{x}}{\mathrm{argmin}}\left({\displaystyle \prod _i{g}_i\left(x_i\right)}\right)$$

(46)

In the integrated system, the factor node constructs the corresponding function to calculate the difference between the predicted measurement and the actual measurement, thereby obtaining the estimate of the state variable and the cost function, that is:

$$f_i\left(x_i\right)=L\left[g_i\left(x_i\right)-z_i\right]$$

(47)

where $z_i$ is the actual measurement value obtained by the sensor, and $L\left(\cdot \right)$ is the cost function.

FGs model the optimal estimation problem in a graphical form and solve for state estimates based on the MAPE criterion. During the optimization process, as the system operates over time, the scale of the graph inevitably increases, leading to a decline in the real-time performance of graph optimization. Therefore, a sliding window mechanism is introduced to balance accuracy and efficiency, thereby enhancing the system real-time capabilities. Table 5 summarizes the performance characteristics of the FG method.

The FG algorithm has been widely applied in various domains, including single GNSS positioning, GNSS and INS integrated positioning, ambiguity resolution, and robust estimation (Watson et al., 2018; Wen et al., 2021; Gao et al., 2022; Zhang et al., 2024). The FG method has been optimized and improved for urban environments to enhance positioning accuracy (Lyu et al., 2021; Bai et al., 2022; Liao et al., 2023; Zhang et al., 2023), and it has been demonstrated that under certain conditions, the FG algorithm achieves higher solution accuracy and robustness compared to the EKF.

2.5. Artificial Intelligence (AI) Method

In addition to the aforementioned fusion methods, machine learning and deep learning have found typical applications in practice, demonstrating excellent performance in fields such as pattern recognition and image processing. Their learning and training concepts have also been applied to position estimation in multi-source integrated navigation, particularly in scenarios where explicit mathematical models cannot be established between navigation source observations and system states. Such as in WLAN positioning, carrier positions cannot be directly calculated from signal strength, and in visual cooperative target localization, carrier positions cannot be directly derived from acquired image tags. However, through offline learning of signal strength or image tags, a learning model can be established between system states and observational information. Thus, in real-time positioning, carrier positions can be directly computed via image tags. By inputting acquired observational data into the trained learning model, system states can be output, enabling multi-source information fusion processing. Table 6 summarizes the advantages and disadvantages of AI methods.

The AI-based algorithms in multi-source data fusion mainly include Neural Networks (NN), Support Vector Machines (SVM), Hidden Markov Models (HMM), Decision Trees, etc. (McClelland et al., 1987; Vapnik, 1999; Zhou, 2016).

NN is one of the most classical AI algorithms. They can be classified based on network topology and information flow direction. When categorized by topological structure, NN models can be divided into hierarchical structures, interconnected structures, and sparsely connected structures according to the connection methods between neurons. Based on the direction of internal information flow, they can be classified as feedforward or feedback neural networks. The learning methods of artificial NN mainly fall into three categories: supervised learning, unsupervised learning, and rote learning. Among them, the classic Backpropagation (BP) neural network belongs to the purely hierarchical type with a feedforward information flow direction and employs supervised learning. Its typical three-layer network model is shown in Figure 7.

The vectors for each layer are denoted as X, Y, and O, respectively. In a multi-source fusion navigation system, the data vector $X=\left[\begin{array}{ccccc}{x}_0& x_1& \dots & x_{n-1}& x_n\end{array}\right]$ refers to the observation information from different sources. The hidden layer vector is denoted as $Y=\left[\begin{array}{ccccc}{y}_0& y_1& \dots & y_{n-1}& y_n\end{array}\right]$ , and the output layer vectors are denoted as $O=\left[\begin{array}{ccccc}{o}_1& \cdots & o_k& \cdots & o_l\end{array}\right]$ . The output can represent the system state, such as position, or large-scale location identifiers. The connection observation noise and weight matrix are denoted as $V=\left[\begin{array}{ccccc}{v}_1& \cdots & v_k& \cdots & v_m\end{array}\right]$ and $W=\left[\begin{array}{ccccc}{w}_1& \cdots & w_k& \cdots & w_l\end{array}\right]$, respectively.

BP neural network algorithm consists of a training process and a learning process, as shown in Figure 8。

Scholars have conducted research on AI-based adaptive algorithms, including fuzzy control adaptive algorithms and neural network adaptive algorithms (Bian et al., 2004). For instance, to address the impact of random disturbances in underwater environments on systems, a Radial Basis Function (RBF) neural network was employed to assist FF for information fusion (Li et al., 2011). In response to the challenges of high costs and significant susceptibility to meteorological conditions in existing high-precision satellite navigation systems for agricultural machinery in farmlands, Yu et al. (2021) proposed a D-S-CNN-based multi-sensor fusion automatic navigation method for farmland applications. These existing AI algorithms require a large amount of training data and pre-scaled training of the system, which makes it difficult to ensure computational capacity and real-time performance. Therefore, they are mostly used for post-computation.

2.6. Methods Based on Uncertain Reasoning

Some fusion methods are not suitable for establishing a mathematical model between measurements and system states, nor can they directly estimate the state. However, they can be used to assess the reliability of navigation sources and determine large-scale locations. For example, the identification of specific large-scale locations such as classrooms, cafes, conference rooms, offices, and residences can be achieved by combining observations from GNSS, WLAN sources, and inertial sources to perform uncertainty reasoning and obtain large-scale location information. Therefore, a multi-source information fusion method based on uncertainty reasoning is proposed. This method requires extracting evidence from measurement data and using relevant knowledge (mainly expert knowledge) to gradually derive conclusions from the evidence or to verify the credibility of specific information. The currently mainstream uncertainty reasoning methods mainly include subjective Bayesian estimation, evidence reasoning, and fuzzy reasoning (Zhang, 1996). This method can also adjust the weights of each measurement data in the final fusion result to achieve precise identification of locations.

Uncertainty reasoning includes symbolic reasoning and numerical reasoning. The former, such as Endorsement Theory, is characterized by minimal information loss during the reasoning process but involves a large computational burden. The latter, such as Bayesian reasoning and Dempster-Shafer theory of evidence, is characterized by ease of implementation but involves some degree of information loss during the reasoning process. As uncertainty reasoning methods are fundamental tools for target recognition and attribute fusion, the subjective Bayes method and the Dempster-Shafer (D-S) evidence theory are two commonly used uncertainty reasoning approaches. For specific applications, please refer to related literature. To facilitate future use, Table 7 summarizes the advantages and disadvantages of uncertainty reasoning methods.

3. Representative Applications of Partial Fusion Methodology in Navigation and Positioning

As the author primarily focuses on research in the field of navigation and location services, only typical applications of some fusion methods in this domain are listed. Contemporary research has demonstrated significant advancements in rapid high-precision navigation for complex urban environments through sensor fusion strategies. Scholars have systematically explored three principal approaches: multi-sensor integration, functional model optimization, and stochastic modeling enhancement. Noteworthy implementations include:

(1) Multi-constellation Hybridization: Wang et al. (2022) developed a train-integrated navigation system leveraging full-state fusion of multi-constellation GNSS and INS. This architecture demonstrates enhanced robustness against signal degradation in metropolitan corridors.

(2) Advanced Tight Coupling Methods: Zhu et al. (2023) proposed a MEMS IMU error mitigation framework combining GNSS RTK carrier phase, TDCP observations, and INS data through adaptive KF. Li et al. (2017) augmented this paradigm with carrier motion constraints, significantly improving solution availability during prolonged GNSS outages.

(3) Precision Positioning Enhancements: Gao (2016) established that multi-GNSS PPP/INS tight coupling reduces PPP convergence time by 40-60% while maintaining centimeter-level accuracy. Urban navigation breakthroughs by Li et al. (2018) revealed that single-frequency RTK/INS integration achieves dual-frequency RTK performance (2-3 cm horizontal RMS) under typical urban multipath conditions.

(4) Multi-Sensor Constraint Integration: Liu (2021) incorporated odometric measurements and zero-velocity constraints into GNSS/INS fusion, demonstrating 75% improvement in INS drift suppression during 60-second GNSS outages. Indoor-outdoor continuity solutions by Yu et al. (2020), Cao (2021), and Yuan et al. (2023) achieved seamless decimeter-level positioning through UWB/GNSS/INS tight coupling.

(5) Resilient Navigation Architectures: Li et al. (2023) engineered a vision-INS-UWB hybrid system where UWB supplants GNSS in denied environments, reducing visual-inertial odometry drift by 32% in extended indoor operations. Emerging multi-sensor SLAM frameworks (Wang et al., 2022) integrate lidar, vision, and inertial data with crowd-sourced mapping, showing particular promise for dynamic urban canyon navigation.

This methodological evolution highlights the critical role of partial fusion techniques in balancing computational complexity with positioning integrity across heterogeneous operational environments.

4. Summary of Features and Application Scenarios for Multi-Source Fusion Methods

The introduced methods for multi-source data fusion processing show differences in terms of optimization criteria, fundamental principles, mathematical models, prior information, number of observations, and application scenarios. To enable users to make targeted selections of different fusion methods when developing integrated systems, we have summarized the main characteristics and applicable scenarios of the fusion methods introduced, as shown in Table 8.

5. Conclusions

This paper provides a detailed overview of various algorithms corresponding to multi-source fusion processing methods. It summarizes the fundamental principles of these algorithms and briefly introduces their mathematical models, key characteristics, and application scenarios, offering significant theoretical and technical support for intelligent navigation, unmanned driving, autonomous navigation, and related fields. Due to limitations in theoretical understanding and technical conditions, all existing fusion algorithms exhibit shortcomings. Currently, no single fusion algorithm can fully meet the requirements of multi-source integrated navigation systems. Therefore, appropriate fusion algorithms must be selected based on practical needs and application contexts. The historical development of these fusion algorithms reveals their interdisciplinary nature, combining theories and methodologies from integrated navigation, GNSS data processing, satellite geodesy, probability theory and mathematical statistics, computer science, statistics, and artificial intelligence. Consequently, multi-source integrated navigation algorithms should not be confined to traditional positioning and navigation approaches. Instead, they should continuously incorporate insights from other disciplines, foster mutual learning and advancement across fields, and generate innovative theories and methods through interdisciplinary integration. This evolution aims to deliver high-precision, high-reliability positioning, navigation, and timing services across all temporal and spatial domains - representing the future development trend of multi-source integrated navigation systems.