**1. Introduction**

The current operating Global Navigation Satellite System (GNSS) includes the Global Positioning System (GPS), Globalnaya Navigazionnaya Sputnikovaya Sistema (GLONASS), BeiDou (BDS), and Galileo. The precise position, velocity, and time (PVT) generated from real-time kinematics (RTK) are crucial for engineering surveys, fleet monitoring, intelligent transportation systems, geographical information systems, guidance and control, etc. [1,2]. By 2024, more than 110 satellites with different frequencies are expected to be accessible for multi-GNSS. Compared with a single constellation, multiple constellations and frequencies could improve the accuracy, continuity, availability, and integrity significantly [3], while enhancing the robustness against noise and outliers [4].

Multi-GNSS applications usually apply loosely and tightly coupled models for data curation [5,6]. Both of them could achieve similar performance [7] in precise point positioning [8] and multi-sensor fusion such as integrating GNSS with the Inertial Navigation

**Citation:** Liu, J.; Liu, T.; Ji, Y.; Sun, M.; Lyu, M.; Xu, B.; Lu, Z.; Xu, G. A Robust Nonlinear Filter Strategy Based on Maximum Correntropy Criterion for Multi-GNSS and Dual-Frequency RTK. *Remote Sens.* **2022**, *14*, 4578. https://doi.org/ 10.3390/rs14184578

Academic Editor: Giuseppe Casula

Received: 30 July 2022 Accepted: 9 September 2022 Published: 13 September 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

System, Simultaneous Location and Mapping, Lidar, and 5G [9–15]. Although the tightly coupled model shares one common pivot satellite for each constellation, it is challenging to improve the model strength due to the presence of double difference inter-system bias (DISB) originating from receivers [16]. In this study, the loosely coupled model is used due to its usability [17].

Even if dual-frequency (DF) data is available [18–20], the ambiguity resolution is still impacted in urban environments due to signal blockage and interference [18–21]. The traditional extended Kalman filter is sensitive to outliers and noise [22], as it is based on the second-order statistics (e.g., variance, correlation, etc.) that exist in the Gaussian noise assumption and minimum mean square error (MSE) criterion [23]. Thus, the following limitations still need to be addressed: (1) The unmodeled non-Gaussian noise and heavytail noise originated not only from the outliers and gross errors but also the missed data due to signal blockage and deformation [24]. (2) The lack of complete prior knowledge of system dynamics and observations may eventually cause divergence [25]. Extensive attempts have been made to address these problems. One solution is the particle filter (PF), which is capable of estimating the posterior probability density function (PDF) by massive particles [26]. Multi-model filters such as the Gaussian sum filter (GSF) are another solution, which parallelly implement and interactively combine several filters to estimate the system states. A computationally economical filter was studied by integrating the robust M estimation into the GSF framework [27]. Heavy-tailed distribution-based filters and the *H*∞ filter have also been studied in [28]. However, attention still needs to be paid, as most of them are deficient in universality and efficiency [29].

Recently, correntropy has received growing attention in signal processing, posture estimation, and machine learning [30–33]. It is a measurable metric of local similarity and is established on Gaussian kernel functions. Specifically, higher dimensional data matching and error detection can be realized, as Gaussian kernel functions enrich data features by transforming the observation into the Hibert space with higher dimensions. Thus, the maximum correntropy between the input and output is robust against various types of noise and even arbitrarily large outliers, and can be achieved according to the maximum correntropy criterion (MCC). Although the kernel bandwidth (KBW) is the key parameter for implementing the MCC, existing studies mostly treat it as an exogenous parameter from empirical experiments, rather than an endogenous variable of the system [33–36]. In this paper, a nonlinear adaptive Kalman filter (KF) based on MCC with adaptive KBW is proposed. The Gaussian hypothesis and MSE criterion are further relaxed, which aims to improve the adaptability and robustness [37]. The main contributions of the proposed method are: (1) The AMC-KF is proposed as a new robust nonlinear filter to improve the precision and robustness of multi-GNSS DF RTK [38]. (2) The DF data-aided AR method is proposed to fix ambiguities with the wide-lane and ionosphere-free combinations. The wide-lane pseudo-range is introduced for medium and long baseline RTK. (3) A nonlinear filter strategy is designed by integrating the DF data-aided AR method into the proposed AMC-KF. The test results show the significant superiority of the proposed strategy with various baselines.

The remainder of this paper is organized as follows: In Section 2, the double difference (DD) RTK model is introduced, followed by the derivation of the proposed DF data-aided AR method. The loosely coupled model for multi-GNSS is also presented. In Section 3, correntropy is introduced and the AMC-KF is elaborated in detail. Moreover, the derivation of the adaptive KBW is outlined. In Section 4, the performance of the proposed filter strategy is demonstrated. Datasets collected from the continuously operating reference stations (CORS) network of Australia are used for the short baseline test, adaptive strategy test, and the long baseline test. Finally, some conclusions are given in Section 5.

#### **2. RTK Mathematic Model**

⎧ ⎪⎪⎪⎪⎪⎪⎨

⎪⎪⎪⎪⎪⎪⎩

*2.1. Constrained Loosely Coupled Model*

The DD carrier phase measurement *ϕ* (in cycles) with wavelength *λ* and code measurement *P* (in meters) is defined as follows [1]:

$$\begin{cases} \lambda \nabla \Delta \rho = \nabla \Delta \rho - \lambda \nabla \Delta N + \nabla \Delta T - \nabla \Delta I + \nabla \Delta \varepsilon\\ \nabla \Delta P = \nabla \Delta \rho + \nabla \Delta T + \nabla \Delta I + \nabla \Delta \xi \end{cases} \tag{1}$$

where ∇Δ is the DD operator; *ρ*, *T*, *I,* and *N* are the receiver-satellite geometric range, tropospheric delay, ionospheric delay, and carrier phase ambiguity, respectively; and *ξ* and *ε* are the unmodeled errors including multipath noise, system noise, etc. The covariance matrix of the system state, system noise, measurement noise, and design matrices can be defined as *P*, *Q, R,* and *H*. The filter-state vector consists of positioning information and the DD ambiguities can be expressed as [7]:

$$\mathbf{x} = \begin{bmatrix} \mathbf{x}\_{\boldsymbol{\upmu}}^{G} & \mathbf{x}\_{\boldsymbol{\upmu}}^{G} & \mathbf{x}\_{\boldsymbol{\upmu}}^{G} & \mathbf{x}\_{\boldsymbol{\upmu}}^{G} & \mathbf{x}\_{\boldsymbol{\upmu}}^{G} & \mathbf{x}\_{\boldsymbol{\upmu}}^{G} & \nabla \Delta \mathbf{N}^{G} & \nabla \Delta \mathbf{N}^{E} \end{bmatrix}^{T}$$

where the superscripts '*G*' and '*E*' represent GPS and Galileo, and the subscripts *e*, *n*, and *u* represent different directions. The *diag*(·) is a diagonal matrix. The KF recursive process is defined as [1]:

$$\begin{aligned} \hat{\mathbf{x}}\_{k|k-1} &= \mathbf{F}\_{k,k-1} \hat{\mathbf{x}}\_{k-1} \\ \mathbf{P}\_{k|k-1} &= \left( \mathbf{F}\_{k,k-1} \mathbf{P}\_{k-1} \mathbf{F}\_{k,k-1}^{\mathrm{T}} + \mathbf{F}\_{k-1} \mathbf{Q}\_{k-1} \mathbf{F}\_{k-1}^{\mathrm{T}} \right) \text{prediction} \\ \mathbf{P}\_{k} &= (\mathbf{I} - \mathbf{K}\_{k} \mathbf{H}\_{k}) \mathbf{P}\_{k|k-1} \\ \mathbf{K}\_{k} &= \mathbf{P}\_{k|k-1} \mathbf{H}\_{k}^{\mathrm{T}} (\mathbf{H}\_{k} \mathbf{P}\_{k|k-1} \mathbf{H}\_{k}^{\mathrm{T}} + \mathbf{R}\_{k})^{-1} \\ \hat{\mathbf{x}}\_{k} &= \hat{\mathbf{x}}\_{k|k-1} + \mathbf{K}\_{k} (\mathbf{z}\_{k} - \mathbf{H}\_{k} \hat{\mathbf{x}}\_{k|k-1}) \end{aligned} \tag{2}$$

where *K*, *F* are the gain and state transform matrix; and *z<sup>k</sup>* and *x<sup>k</sup>* are the measurements and state that needs to be estimated at the *k*th epoch. As *x<sup>G</sup>* = *x<sup>E</sup>* once successfully located, the following constraint can be attached [39]:

$$\underbrace{\begin{bmatrix} \mathbf{I}\_{3\times3} & -\mathbf{I}\_{3\times3} \end{bmatrix}}\_{D} \underbrace{\begin{bmatrix} \mathbf{x}\_{G} \\ \mathbf{x}\_{E} \end{bmatrix}}\_{X} = \underbrace{\begin{bmatrix} \mathbf{0}\_{3\times3} \\ \mathbf{0}\_{M} \end{bmatrix}}\_{M} \tag{3}$$

where *I* represents the unit matrix, and *D* and *M* are the constraint matrices. The prediction step in (2) can be developed to constrain KF as follows:

$$\begin{cases} \mathbf{x}\_{k|k-1} = \mathbf{x}\_{k|k-1} - \mathbf{D}\_{k}^{T} \left(\mathbf{D}\_{k} \mathbf{D}\_{k}^{T}\right)^{-1} \left(\mathbf{D}\_{k} \mathbf{D}\_{k}^{T}\right)^{-1} \left(\mathbf{D}\_{k} \mathbf{x}\_{k|k-1} - \mathbf{M}\_{k}\right) \\\ \mathbf{P}\_{k|k-1} = \quad \left(\mathbf{I} - \mathbf{D}\_{k}^{T} \left(\mathbf{D}\_{k} \mathbf{D}\_{k}^{T}\right)^{-1} \mathbf{D}\_{k}\right)^{T} \mathbf{P}\_{k|k-1} \left(\mathbf{I} - \mathbf{D}\_{k}^{T} \left(\mathbf{D}\_{k} \mathbf{D}\_{k}^{T}\right)^{-1} \mathbf{D}\_{k}\right) \end{cases} \tag{4}$$

#### *2.2. DF Data-Aided AR*

Defining *N*<sup>1</sup> and *N*<sup>2</sup> as the carrier ambiguities on frequencies *f* <sup>1</sup> and *f* 2, respectively, the *ρ* on different frequencies can be formed by ∇Δ*ϕ* as follows:

$$\begin{cases} \begin{array}{l} \nabla \Delta \rho\_1 = (\nabla \Delta \varrho\_1 + \nabla \Delta N\_1) \lambda\_1 - \frac{A}{f\_1^2} \\ \nabla \Delta \rho\_2 = (\nabla \Delta \varrho\_2 + \nabla \Delta N\_2) \lambda\_2 - \frac{A}{f\_2^2} \end{array} \tag{5}$$

where *Ne* is the number of electrons in unit area *A* = 40.3, *<sup>S</sup> Neds*. The wide-lane ambiguity ∇Δ*Nw* with wavelength *λ<sup>w</sup>* can be expressed as:

$$
\nabla \Delta N\_w = \nabla \Delta q\_1 - \nabla \Delta q\_2 - \frac{1}{\lambda\_w} \left( \nabla \Delta \rho - \nabla \Delta T - \frac{f\_1}{f\_2} \nabla \Delta I - \nabla \Delta \varepsilon \right) \tag{6}
$$

For the atmosphere errors in the above equation, ∇Δ*I* cannot be ignored for long baselines due to the spatial difference. ∇Δ*T* varies from 2 to 20 m depending on the satellite elevation, but a higher cutoff angle reduces data utilization [40]. Thus, the timeliness of the fixed solution could not be guaranteed, which is crucial for dynamic RTK once the satellites are available [41–44]. The proposed DF data-aided AR method is summarized as follows:

$$\begin{cases} \nabla \Delta \mathbf{N}\_{\mathcal{w}} = (\nabla \Delta q\_1 - \nabla \Delta q\_2) - \frac{(f\_1 \nabla \Delta P\_1 + f\_2 \nabla \Delta P\_2)}{\lambda\_{\mathcal{w}} (f\_1 + f\_2)} + \nabla \Delta \varepsilon\\ \nabla \Delta \mathbf{N}\_1 = \frac{1}{m \lambda\_1 - n \lambda\_2} [\nabla \Delta q \rho\_{IF} - m \nabla \Delta q\_1 + n \nabla \Delta q\_2 - n \lambda\_2 \nabla \Delta \mathbf{N}\_{\mathcal{w}}] \\\ \nabla \Delta \mathbf{N}\_2 = \nabla \Delta \mathbf{N}\_1 - \nabla \Delta \mathbf{N}\_{\mathcal{w}} \end{cases} \tag{7}$$

where <sup>∇</sup>Δ*ϕ*IF is the ionosphere-free combination for carrier observations, *<sup>m</sup>* <sup>=</sup> *<sup>f</sup>* <sup>2</sup> 1 *f* 2 <sup>1</sup> <sup>−</sup>*<sup>f</sup>* <sup>2</sup> 2 and *<sup>n</sup>* <sup>=</sup> *<sup>f</sup>* <sup>2</sup> 2 *f* 2 <sup>1</sup> <sup>−</sup>*<sup>f</sup>* <sup>2</sup> 2 , and ∇Δ*ε* represents the noise. The above method introduces the code wide-lane combination to invert ∇Δ*N*<sup>1</sup> and the derivation can be found in Appendices A and B. The proposed method shows the following merits: (1) The positioning accuracy and AR success rate are improved as the influence of ∇Δ*I* error is eliminated. (2) The method is applicable for medium and long baselines as the limitation of geometry distance *ρ* is eliminated. (3) The ∇Δ*N* on each frequency can be inverted directly from the high precise ionosphere-free and wide-lane measurements. Moreover, the following moving average with *n* epochs is also adopted to reduce the effect result from noise in ∇Δ*P* [45]:

$$\nabla \Delta \hat{\mathcal{N}}\_w(n) = \frac{\sum\_{i=k}^{k+n} Z^i \nabla \Delta \mathcal{N}\_w^i}{\sum\_{\substack{k+n\\i=k}}^{k+n} Z^i} \tag{8}$$

where *k* is the start epoch of the observation arc without cycle slip. *Zi* is the weight of *<sup>i</sup>*th epoch. The <sup>∇</sup>Δ*N*<sup>ˆ</sup> *<sup>w</sup>*(*n*) could be fixed to *round* <sup>∇</sup>Δ*N*<sup>ˆ</sup> *<sup>w</sup>*(*n*) by the integer rounding method [46] if the differential residual of <sup>∇</sup>Δ*N*<sup>ˆ</sup> *<sup>w</sup>* between adjacent epochs meets the following constraint:

$$\left| \nabla \Delta \hat{\aleph}\_w(n) - \nabla \Delta \hat{\aleph}\_w(n-1) \right| < 0.1 \text{ cycles}$$

Then, the corresponding ∇Δ*N*1, ∇Δ*N*<sup>2</sup> and precise carrier measurements could be formed by Equation (7). In the following test, the moving average window width *n* is set to be five epochs.

#### **3. AMC-KF**

The obtained ∇Δ*N*<sup>1</sup> and ∇Δ*N*<sup>2</sup> could be applied to form precise ∇Δ*ϕ* in Equation (1), and used for state estimating in the nonlinear filter.
