A NLOS Ranging Error Mitigation Method for 5G Positioning in Indoor Environments

1. Introduction

Temporal and spatial information is essential to human existence and development. Location-Based Services (LBS) utilize positioning and navigation technologies to provide users with a series of services based on temporal and spatial locations, playing an increasingly important role in various emerging fields such as smart cities, intelligent transportation, and industrial internet. Ubiquitous and trustworthy location services require the support of advanced positioning technologies. In outdoor environments, GNSS can provide high precision location information. However, due to the low penetrability of satellite signals, it is not possible to offer users high precision location services in non-exposed spaces such as indoor and underground environment [1,2]. Currently, positioning technologies in satellite denied environments worldwide include Bluetooth, UWB, WLAN, and mobile communication networks [3]. Among these, mobile communication network is the most widely used wireless networks for ground coverage. Especially 5G networks can provide high-precision location services based on mobile communication networks has become a research hotspot [4,5,6,7,8]. The primary technical bottleneck in improving 5G indoor positioning accuracy is the NLOS effect [9]. In complex environments such as indoors and underground, wireless signal propagation is severely obstructed and reflected, resulting in NLOS phenomena [10]. Since communication services in mobile communication networks are insensitive to delays, NLOS phenomena have a minor impact on communication quality. However, for positioning, NLOS can lead to increased positioning observation errors. The errors caused by NLOS may be much larger than the errors introduced by the ranging capability itself. Therefore, researching NLOS error suppression techniques in mobile communication networks is of great significance for achieving high-precision positioning.

Recently a large number of scholars have conducted extensive research on the NLOS error generation mechanism and compensation technology. The existing NLOS suppression technology is mainly divided into three categories, as shown in Table 1.

The first category is detection, which involves initially identifying the observations of positioning signals. Observations that may be affected by NLOS are either directly removed or minimized in terms of their weight. The advantage of this approach is its simplicity, while its disadvantage is the lack of classification. [11] proposed the consistency checking method to examine the residual range measurement or navigation solution, which is similar to the receiver autonomous integrity monitoring.

The second category is mitigation, which uses multiple features to classify the received positioning signals, and performs error suppression based on the observation feature classification results. Technologies such as machine learning and support vector machines have been widely used to classify and feature extraction of GNSS and wireless communication network observations [14,15,16,17]. [18] proposed the reception of NLOS signals based on Vector tracking loop, and initially proved the effectiveness of the proposed method through field tests. However, there is currently a lack of comprehensive theoretical derivation and performance limit analysis for this type of method.

The third type is correction. In this type of technique, NLOS paths are considered as useful signals. Utilizing known three-dimensional maps for assisted positioning, potential reflection paths are inferred, and error correction is performed accordingly [19,20,21,22]. This is equivalent to obtaining a new observation quantity that is not affected by non-line-of-sight errors. This type of method can ideally remove the impact of NLOS errors on positioning, and can even bring additional improvements in performance. However, the computational complexity is very high, and it is difficult to meet the real-time requirements. In addition, the deployment on hardware platform is also a significant challenge.

The existing NLOS mitigation methods suffer from insufficient feature extraction of wireless signals, lack consideration of residual NLOS measurement errors during the position estimation process, and poor stability in the mixed positioning estimation results [23,24,25]. To address these issues, this paper innovatively proposed a NLOS error identification, suppression, and position estimation method for 5G positioning terminals. The primary contributions of this paper are outlined as follows:

1) A method for discriminating and classifying NLOS observations based on decoupling the features of wireless signals is proposed. This method initially divides the indoor propagation channels into multiple categories and utilizes the results of channel classification to assess the severity of NLOS effects. Subsequently, it establishes selection criteria for filtering the ranging results of 5G positioning signals.

2) A position estimation method based on Square-Root Unscented Kalman Filter (SRUKF) is proposed. The square root of the state covariance is introduced, thereby avoiding the need for covariance reconstruction at each iteration. This approach offers advantages such as numerical stability and positive semi-definiteness of the state covariance. During the fusion process of mixed observations, it ensures the numerical stability of the filter’s position estimation output.

3) We perform a comprehensive comparison of the existing advanced methods and build a 5G localization system to support the validation of the proposed method in real complex indoor environment. Both the simulation and experimental results demonstrate the superiority of the proposed NLOS mitigation and location estimation method over existing methods.

The remainder of this paper is organized as follows. In Section 2 the 5G positioning reference signal measurement model is given and related problems are described. In Section 3 a NLOS mitigation and location estimation method for 5G positioning is proposed which is composed of the derivation of positioning performance limits, measurement identification and location estimation method for NLOS mitigation. In Section 4 the performance of the proposed NLOS error mitigation algorithm is evaluated through the simulations as well as the experiment test. Finally, the discussion of this article are drawn in Section 5.

2. Problem Description and System Model

This section introduces the 5G NR time-frequency resource mapping scheme. The 5G positioning reference signal (PRS) has an abundance of time-frequency resources which is expected to show better performance. The pseudorandom code of the PRS signal is generated by gold sequences and OFDM modulation is employed to achieve reliable data transmission in complex scenarios with a scalable OFDM numerology. The physical signal’s transmissions are organized into signal frames, subframes, and slots in the time domain. Each signal frame has a duration of 10ms and consists of 10 subframes with a subframe duration of 1ms. A subframe is formed by one or multiple slots. The resource element consists of one subcarrier in one OFDM symbol, which is the smallest physical time-frequency resource unit. The physical resource block is defined as 12 continuous subcarriers in one time slot.

In this study, factors such as positioning accuracy and resource occupancy rate are considered comprehensively. The parameter of 30KHz subcarrier spacing was used in our experiments and thus, one subframe has 2 slots, and each slot has 14 OFDM symbols. Comb-4 is adopted for resource mapping. This paper takes into account a wide range of parameters, including resource occupancy rate and placement accuracy. In our experiments, we selected a 30 kHz subcarrier spacing setting, which results in two slots per subframe with 14 OFDM symbols in each slot. Comb-4 is utilized for the mapping of resources.

In this section, we analyze the 5G PRS as an example. After appropriate adjustments, the proposed method can also be applied to the suppression of NLOS errors in other wireless signals. The baseband received signal at the terminal for 5G positioning can be represented as follows:

$$\begin{array}{cc}\hfill r\left(t\right)=& \sum _{m=1,m\ne n}^M\left(A^m{C}^m\left(t-{\tau }^m\right)cos\left(2\pi \left(f_{IF}^m+f_d^m\right)t+{\phi }^m\right)\right)\hfill \\ \hfill & +{\alpha }^{\mathrm{NLOS}}{A}^n{C}^n\left(i-{\tau }^{\mathrm{LOS},n}-{\tau }^{\mathrm{NLOS}}\right)\hfill \\ \hfill & \times cos\left(2\pi \left(f_{IF}^n+f_d^n\right)t+{\phi }^n\right)+\omega \left(t\right)\hfill \end{array}$$

(1)

Where $r\left(t\right)$ represents the 5G signal sampled at time t, m denotes the identifier of the positioning base station in the network, M represents the total number of positioning base stations. A indicates the amplitude value of the positioning signal, C represents the ranging code sequence, ${\tau }^m$denotes the delay of the ranging code for base station, $f_{IF}$ and $\phi$ denote the frequency and initial phase of the local carrier respectively, $f_d$ represents the Doppler frequency shift, n indicates the identifier of the NLOS base station, ${\alpha }^{\mathrm{NLOS}}$ represents the reflection coefficient of the NLOS signal, ${\tau }^{\mathrm{NLOS}}$ denotes the delay caused by the NLOS effect relative to the direct path signal and $\omega \left(t\right)$ represents additive Gaussian white noise.

Achieving high-precision positioning of 5G PRS signals requires continuous signal tracking at the terminal side. The impact of NLOS effects on the code tracking loop is far greater than that on the carrier tracking loop. Assuming the terminal can completely strip the carrier from the received signal, the outputs of the early, prompt, and L correlators can be represented as follows:

$$E=A^m·R\left({\tau }_e^m+\frac{d}{2}\right)+{\alpha }^{\mathrm{NLOS}}{A}^n·R\left({\tau }_e^{\mathrm{NLOS}}+\frac{d}{2}\right)$$

(2)

$$P=A^m·R\left({\tau }_e^m\right)+{\alpha }^{\mathrm{NLOS}}{A}^n·R\left({\tau }_e^{\mathrm{NLOS}}\right)$$

(3)

$$L=A^m·R\left({\tau }_e^m-\frac{d}{2}\right)+{\alpha }^{\mathrm{NLOS}}{A}^n·R\left({\tau }_e^{\mathrm{NLOS}}-\frac{d}{2}\right)$$

(4)

Where R represents the autocorrelation function, d represents the spacing between the early and late correlators, ${\tau }_e^m={\widehat{\tau }}^m-{\tau }^m$ represents the code phase difference of the direct path signal, ${\tau }_e^{NLOS}={\widehat{\tau }}^{NLOS}-\left({\tau }^{LOS,n}+{\tau }^{NLOS}\right)$ is the code phase difference of the NLOS signal and ${\widehat{\tau }}^m$ represent the estimated values of the ranging results. Here, we mainly consider the scenario of 2D positioning, but it can be extended to 3D positioning as needed. The estimated values of ranging results at the terminal can be modeled as:

$${\widehat{d}}_m=\left\{\begin{array}{cc}{d}_m+{\omega }_m,\hfill & \mathrm{for}\mathrm{LOS}\mathrm{case}\hfill \\ d_m+{\omega }_m+{\eta }_m,\hfill & \mathrm{for}\mathrm{NLOS}\mathrm{case}\hfill \end{array},m=1,2,\cdots ,M\right.$$

(5)

$$d_m=\sqrt{{\left(x_m-x\right)}^2+{\left(y_m-y\right)}^2}$$

(6)

$\left(x_m,y_m\right)$ and $(x,y)$ represents the coordinates of 5G base station and the terminal. ${\widehat{d}}_m$ represents the estimated distance from the mth base station to the terminal. ${\omega }_m$ denotes zero-mean Gaussian noise, ${\eta }_m$ is NLOS error term. The ranging results from M base stations are divided into two categories, ${\Phi }_{LOS}$ and ${\Phi }_{NLOS}$ , representing line-of-sight and non-line-of-sight, respectively.

$${\Phi }_{LOS}=\left\{\left(x_m,y_m,{\widehat{d}}_m\right):{\widehat{l}}_m=1\right\},{\Phi }_{NLOS}=\left\{\left(x_m,y_m,{\widehat{d}}_m\right):{\widehat{l}}_m=-1\right\}$$

(7)

Where ${\widehat{l}}_m$ is the indicator of channel state, 1 represents LOS channel, and -1 represents NLOS channel. In indoor environment the terminal often needs to provide position estimation results in the presence of a mixture of LOS and NLOS ranging results. The optimal set of ranging results available for the user terminal to choose from can be denoted as

$$\Phi ={\Phi }_{LOS}\cup {\Phi }_{NLOS}^{\prime }\mathrm{and}{\Phi }_{NLOS}^{\prime }\subseteq {\Phi }_{NLOS}$$

(8)

Afterwards, there can be various methods to provide position estimation values based on ranging results. However, existing methods have not considered the incomplete identification of NLOS observations or incomplete suppression of NLOS errors. These residual NLOS errors can lead to a significant increase in positioning errors for 5G positioning.

3. The Proposed NLOS Mitigation and SRUKF Based Location Estimation Method for 5G Positioning

In this chapter, a NLOS mitigation and SRUKF based location estimation method (NM-SL) for 5G positioning terminals is proposed. Firstly, the performance limits of TOA estimation in indoor environments with mixed line-of-sight and NLOS channels is derived. Based on this, a strategy for filtering NLOS observations is proposed and a position estimation method for NLOS ranging errors based on SRUKF is introduced. The diagram of the proposed method is illustrated in Figure 1.

3.1. Ziv-Zakai Bound on Time-of-Arrival Estimation with Statistical NLOS Channel Knowledge at the Receiver

In this section the performance limits of TOA estimation in NLOS channels is analyzed. Analytical expressions for the Cramér-Rao Bound (CRB) and Ziv-Zakai Bound (ZZB) of the root mean square error of in different NLOS conditions are derived. These performance limits can serve as effective benchmarks for designing practical TOA estimators. The 5G positioning terminal obtains estimation of TOA, denoted as ${\widehat{\tau }}^m$, from the received signal r(t) in Equation (1). We assume that $\tau$ is unknown and randomly distributed within the sampling time interval, and all NLOS components are observed within this sampling interval. The proposed performance boundary based on ZZB can be derived from the following general equation for root mean square error estimation

$$\mathrm{E}\left\{{\xi }^2\right\}=\frac{1}{2}{\int }_0^{\infty }\phi ·\mathrm{P}\left\{\left|\xi \right|\ge \frac{\phi }{2}\right\}d\phi$$

(9)

where $\xi ={\widehat{\tau }}^m-{\tau }^m$ represents the estimation error of the ranging results in the positioning terminal. The core idea of ZZB is to transform the performance evaluation of an estimation problem into a binary detection problem with equal probability hypotheses.

$$\begin{array}{cc}\hfill & H_1:r\left(t\right)\sim \mathbf{p}\{r\left(t\right)\mid \tau \}\hfill \\ \hfill & H_2:r\left(t\right)\sim \mathbf{p}\{r\left(t\right)\mid \tau +\phi \}\hfill \end{array}$$

(10)

Where $\mathrm{p}\{·\}$ represents following probability distribution, it can be used to obtain the error probability corresponding to the optimal decision rule based on likelihood ratio test, further derived from formula (9) to obtain its performance lower bound. First, estimate the value of ${\tau }^m$, then make a decision between two hypotheses based on the minimum distance criterion, and obtain the performance lower bound by calculating the error probability of the optimal decision rule that minimizes the error probability.

$${\rm Y}\left(r\left(i\right)\right)=\frac{\mathbf{p}\left\{r\right(t)\mid \tau \}}{\mathbf{p}\left\{r\right(t)\mid \tau +\phi \}}\left\{\begin{array}{c}{H}_1:>1\hfill \\ H_2:<1\hfill \end{array}\right.$$

(11)

When ${\tau }^m$ is uniformly distributed within the sampling interval T, then ZZB can be expressed as

$$ZZB=\frac{1}{T}{\int }_0^T\phi (t-\phi )P_{min}\left(\phi \right)d\phi$$

(12)

From the analysis above, it can be seen that the derivation of ZZB for NLOS channels has been transformed into the performance evaluation of a binary detection problem. Deriving ZZB requires evaluating the probability error function $P_{min}\left(\phi \right)$ for decision making between two hypotheses in Equation (11). Under the condition where the channel state information is known, the ideal received signal $\mathrm{r}\left(\mathrm{t}\right)$ can be obtained along with the corresponding channel parameters $\mathit{c}={\left[c_1,c_2,\dots ,c_L\right]}^{\mathrm{T}}$, and the likelihood ratio test results can be computed. This process can be equivalently interpreted as evaluating the corresponding error probability $P_{min}(\phi ,\mathit{c})$ for a binary discrimination system.

$$\begin{array}{cc}\hfill & H_1:r(t-\tau ,\mathit{c})\hfill \\ \hfill & H_2:r(t-\tau -\phi ,\mathit{c})\hfill \end{array}$$

(13)

Where $\mathit{c}$ is known, the error probability can be expressed as

$$P_{min}(\phi ,\mathit{c})=Q\left(\sqrt{{\kappa }_{\mathrm{snr}}\left({\rho }_r(0,\mathit{c})-{\rho }_r(\phi ,\mathit{c})\right)}\right)$$

(14)

Where ${\kappa }_{\mathrm{snr}}=E_{\mathrm{p}}/N_0$ denotes the signal-to-noise ratio, $E_{\mathrm{p}}$ is the average received energy, Q denotes the Gaussian Q-function. The normalized autocorrelation function is defined as follows:

$$\begin{array}{ccc}\hfill {\rho }_r(\phi ,\mathit{c})& =\frac{1}{E_{\mathrm{p}}}{\int }_{-\infty }^{\infty }r(t-\tau ,\mathit{c})r(t-\tau -\phi ,\mathit{c})dt\hfill & \hfill ={\mathit{c}}^T\mathbf{R}\left(\phi \right)\mathit{c}\end{array}$$

(15)

where

$$\mathbf{R}\left(\phi \right)\triangleq \left(\left[\begin{array}{cc}cos\left(2\pi f_c\phi \right)\hfill & sin\left(2\pi f_c\phi \right)\hfill \\ -sin\left(2\pi f_c\phi \right)\hfill & cos\left(2\pi f_c\phi \right)\hfill \end{array}\right]\otimes \mathbf{R}\left(\phi \right)\right)$$

(16)

By replacing $P_{min}\left(\phi \right)$ with $P_{min}(\phi ,\mathit{c})$ in Equation (12), we can obtain the ZZB under the condition of channel state $\mathit{c}$.

$${\left.ZZB\right|}_c=\frac{1}{T}{\int }_0^T\phi (T-\phi )P_{min}(\phi ,\mathit{c})d\phi$$

(17)

On the basis of (17), by further deriving the idea of averaging on the vector c, the unconditionally restricted ZZB can be obtained

$$ZZB={\left.\int f_{\mathit{c}}\left(\mathit{c}\right)ZZB\right|}_{\mathit{c}}d\mathit{c}=\frac{1}{T}{\int }_0^T\phi (T-\phi ){\overline{P}}_{min}\left(\phi \right)d\phi$$

(18)

The average error probability function can be expressed as

$${\overline{P}}_{\min }\left(\phi \right)=\int f_{\mathit{c}}\left(\mathit{c}\right)P_{min}(z,\mathit{c})d\mathit{c}$$

(19)

It should be noted that the premise of obtaining Equation (18) is that the positioning terminal has accurate channel state information. Therefore, when the positioning terminal receives signals completely synchronized, ZZB can serve as a lower bound for the root mean square error of TOA estimation. Furthermore, Equation (19) represents the average bit error probability using receiver coherent detection with perfect channel state information in the presence of NLOS and multipath. Many scholars have studied the coupling relationship between BEP and ZZB under different conditions [26]. Here, we adopt a method similar to [27] to solve for ${\overline{P}}_{min}\left(\phi \right)$

$$Q\left(x\right)=\frac{1}{\pi }{\int }_0^{\pi /2}exp\left\{-\frac{x^2}{2{sin}^2\left(\theta \right)}\right\}d\theta$$

(20)

Substituting Equation (14) and Equation (20) into Equation (19) can obtain the expression for the error probability function as follows

$$\begin{array}{cc}\hfill {\overline{P}}_{min}\left(\phi \right)& ={\mathrm{E}}_{\mathit{c}}\left\{P_{min}(\phi ,\mathit{c})\right\}\hfill \\ \hfill & =\frac{1}{\pi }{\int }_0^{\pi /2}{E}_{\mathit{c}}\left\{exp\left\{-\frac{{\kappa }_{\mathrm{snr}}{\mathit{c}}^T\left[{\mathbf{R}}_2\left(0\right)-{\mathbf{R}}_2\left(\phi \right)\right]\mathit{c}}{2{sin}^2\left(\theta \right)}\right\}\right\}d\theta \hfill \\ \hfill & =\frac{1}{\pi }{\int }_0^{\pi /2}{\Phi }_{\rho }\left(-\frac{{\kappa }_{\mathrm{snr}}}{2{sin}^2\left(\theta \right)}\right)d\theta \hfill \end{array}$$

(21)

Where ${\Phi }_{\rho }$ represents the moment generating function (MGF). The specific derivation process can be referred to [28]. Here, we only need to evaluate an integral with a finite limit. An effective method for evaluating the integral in the form of Equation (21) can be found in [29]. Therefore, under the condition of known accurate channel state information in an Additive White Gaussian Noise (AWGN) channel, the performance bound ZZB for wireless localization can be expressed as:

$$ZZB=\frac{1}{T}{\int }_0^T\phi (T-\phi )Q\left(\sqrt{{\kappa }_{\mathrm{snr}}\left(1-{\rho }_{\mathrm{p}}\left(\phi \right)\right)}\right)d\phi$$

(22)

3.2. 5G NLOS Measurement Identification and Classification Method

When the available LOS signals are not sufficient to support localization, it becomes necessary to design a criterion for filtering ranging results to ensure that only beneficial NLOS signals are included in the positioning computation. The NLOS error at the 5G positioning terminal can be represented as

$${\mathrm{Error}}_{NLOS}=\frac{10\left(\sqrt{{\vartheta }_{nlos}}-1\right)}{{ϵ}_{nlos}}{log}_{10}\left(\frac{{max}_{1<i\le L}\left(a_i^2{\tau }_i^n\right)}{a_{DP}^2{\tau }_{DP}^n}\right)$$

(23)

Where ${\vartheta }_{\mathrm{nlos}}$ represents the dielectric constant associated with obstacles in the propagation channel, ${ϵ}_{\mathrm{nlos}}$ represents the attenuation constant associated with obstructing objects in the channel. L represents the number of paths for multipath, $a_i$ and ${\tau }_i$ respectively represent the amplitude and delay, the subscript DP stands for the direct path signal. It can be seen that the errors introduced by NLOS effects are closely related to the channel conditions. The key to select NLOS observations lies in extracting typical classification features from the waveform of the received signal. The extracted features must be strongly correlated with changes in the NLOS channel state.

To ensure the ease of implementation of the proposed method, we selected six common features from the 5G signal: Total Energy (TE), Maximum Peak Value (MPV), Signal-to-Noise Ratio (SNR), Rise Time (RT), Root Mean Square Delay Spread (RDS), and Kurtosis (KUR). Their definitions are shown in equations (24) - (29) respectively.

$$\zeta =\sum _{i=1}^N{\left|r\left(t_i\right)\right|}^2$$

(24)

$$A_{max}=max\left\{\left|r\left(t_i\right)\right|\right\}$$

(25)

$${\kappa }_{snr}=10{log}_{10}\left(\frac{A_{max}^2}{2{\sigma }_n^2}\right)$$

(26)

$$\begin{array}{c}{\tau }_{\mathrm{rise}}={\tau }_{\mathrm{stop}}-{\tau }_{\mathrm{start}}\\ \left\{\begin{array}{c}{\tau }_{\mathrm{start}}=min\left\{t_i:\left|r\left(t_i\right)\right|\ge 0.1A_{\max }\right\}\hfill \\ {\tau }_{\mathrm{stop}}=min\left\{t_i:\left|r\left(t_i\right)\right|\ge 0.9A_{\max }\right\}\hfill \end{array}\right.\end{array}$$

(27)

$$\lambda =\frac{1}{N{\sigma }_r^4}\sum _{i=1}^N{\left(\left|r\left(t_i\right)\right|-{\mu }_r\right)}^4$$

(28)

$${\tau }_{rms}=\frac{1}{\zeta }\sum _{i=1}^N\left[{\left(t_i-\frac{1}{\zeta }\sum _{i=1}^N\left(t_i{\left|r\left(t_i\right)\right|}^2\right)\right)}^2{\left|r\left(t_i\right)\right|}^2\right]$$

(29)

Where $\mathrm{N}$ is the total number of sampling points of the signal waveform received by the terminal, and i represents the sampling point index. ${\mu }_r=\frac{1}{N}{\sum }_{i=1}^N\left|r\left(t_i\right)\right|,{\sigma }_r^2=$$\frac{1}{N}{\sum }_{i=1}^N{\left(\left|r\left(t_i\right)\right|-{\mu }_r\right)}^2$. When a wireless signal encounters obstacles, its total signal energy will inevitably experience significant attenuation. Additionally, LOS channels obstructed by obstacles can lead to significant fluctuations in RT, RDS, and KUR values. Moreover, areas with severe multipath effects often accompany severe NLOS errors, causing RDS to increase with the increase in multipath classification. Finally, KUR can also serve as an important indicator for assisting in judging the quality of NLOS observations. KUR values decrease significantly with the increase in severity of multipath interference.

The Fuzzy Comprehensive Evaluation is adopted to classify channel states. Let $\widehat{F}={\left\{f_x\right\}}_{x=1}^{\alpha }$ represent the feature set extracted from the signal waveform, and the corresponding distance estimation result be $\widehat{d}$, where $\alpha$ denotes the total number of classification features. The specific method for classifying NLOS channels is as follows:

Find the data closest to real-time measurements in the local sampled dataset. Let the local dataset under channel p be denoted as $D.\widehat{F}$ and $\mathrm{D}$ can be respectively considered as the factor domain and evaluation domain. Assuming the weight values $w_k$ for the above features $f_k$ regarding channel classification, the combination importance of $\alpha$ classification features can be viewed as a fuzzy set $\underset{\sim }{W}$ in the feature domain $\widehat{F}$. The evaluation result is defined as a fuzzy set $\underset{\sim }{S}={\left\{s_x^i\right\}}_{x=1}\left(s_x^i\in [0,1]\right)$ in domain $\mathrm{D}$, which can be obtained through the following fuzzy transformation:

$$\underset{\sim }{S}=\underset{\sim }{W}·\underset{\sim }{R}$$

(30)

$$W={\left\{w_k\right\}}_{k=1}^{\alpha },\mathrm{subject}\mathrm{to}\left\{\begin{array}{c}0\le w_k\le 1\hfill \\ \sum _{x=1}^{\alpha }{w}_k=1.\hfill \end{array}\right.$$

(31)

$\underset{\sim }{R}$ can be expressed in matrix form as $\left[r_{kx}\right]$. In this article, the fuzzy operation operation can be regarded as the inner product of two vectors.

$$s_x^i=\sum _{k=1}^p\left(w_k{r}_{kx}\right)$$

(32)

The sampling data that achieves the maximum value in the evaluation result $\mathrm{S}$ is considered as the data closest to real time measurement. The corresponding $\widehat{F}$ can be represented as $F_{\mathrm{best}}$. A new local dataset $D_{\mathrm{best}}={\left\{F_{\mathrm{best}}^i,\Delta d_{\mathrm{best}}^i,l_{\mathrm{best}}^i,s_{\mathrm{best}}^i\right\}}_{i=1}^p$ can be generated, where $\Delta d^i$ represents the ranging error, $l_{\mathrm{best}}^i$ represents the current true channel condition. We identify the propagation channel of online measurements as one of the sample channels in $D_{\mathrm{best}}$, or a mixture of multiple sample channels. When NLOS is determined, we can further predict the NLOS ranging error using the same method. There is an upper limit on the representative number of channels, the threshold can be expressed as

$$s_{\mathrm{best}}^i<\beta ,\forall i,i\in \{1,2,\dots ,p\}$$

(33)

According to Equation (33) if the volume of real-time test data exceeds the discrimination capability of this method, it is necessary to discard all data to avoid outliers that deviate significantly from the local data. Subsequently, the implementation of the NLOS identification and suppression method is as follows:

$$\left\{\begin{array}{cc}\mathrm{nlos},\tilde{d}=\widehat{d}-\Delta d_{\mathrm{best}}^i,\hfill & \mathrm{when}\widehat{l}=l_{\mathrm{best}}^i<0\hfill \\ los,\tilde{d}=\widehat{d},\hfill & \mathrm{when}\widehat{l}=l_{\mathrm{best}}^i>0.\hfill \end{array}\right.$$

(34)

If the above conditions cannot be met, the positioning channel at this time can be divided into a combination of multiple channels, and these channels all meet the conditions

$$s_{\mathrm{best}}^i\ge \beta ,i\in \{1,2,\dots ,p\}$$

(35)

The NLOS suppression method can be denoted as

$$\left\{\begin{array}{c}NLOS,\tilde{d}=\widehat{d}-\sum \left(\frac{s_{\mathrm{best}}^i}{\Psi }\Delta d_{\mathrm{best}}^i\right),\mathrm{when}\\ \widehat{l}=\sum \left(\frac{s_{\mathrm{best}}^i}{\Psi }{l}_{\mathrm{best}}^i\right)<0\\ \mathrm{LOS},\tilde{d}=\widehat{d},\mathrm{when}\widehat{l}=\sum \left(\frac{s_{\mathrm{best}}^i}{\Psi }{l}_{\mathrm{best}}^i\right)\ge 0,\end{array}\right.$$

(36)

Where

$$\Psi =\sum s_{\mathrm{best}}^i,\mathrm{subject}\mathrm{to}{s}_{\mathrm{best}}^i\ge \beta$$

(37)

3.3. SRUKF Based Location Estimation Method for NLOS Mitigation

In the previous section, we presented a method for NLOS mitigation for each base station’s positioning signal. In this section, we select some NLOS observations based on the previous section for fused position estimation. To improve the numerical stability of the output results of the position estimator, we introduce the square root of the state covariance on the basis of the UKF. In the implementation of SRUKF, the square root of the state covariance is directly propagated, avoiding the need to reconstruct the covariance at each iteration. Additionally, the square root form has the additional advantages of numerical stability and positive semi-definiteness of the state covariance. It is essential to ensure the stability of the filter during the fusion process of LOS and NLOS observations. The implementation of the proposed method is shown in Figure 2.

The initial system state follows a Gaussian distribution with mean m and variance P. From the system state matrix we can get 2L+1 sigma points. [30] presented the way to calculate sigma point. The nonlinear function is used to propagate these sigma vectors.

$${\mathcal{Y}}_i=h\left(x_i\right)i=0,\dots ,$$

(38)

The covariance of the 5G LOS and NLOS mixed positioning state prediction error can be expressed as

$$P^{-}\left[t_k\right]=FP\left[t_k\right]F^T+Q\left[t_k\right]$$

(39)

Where $Q\left[t_k\right]$ is Gaussian process noise matrix at $t_k$ and $P\left[t_k\right]$ denotes the estimate error covariance and its square root form $R\left[t_k\right]$

$$P\left[t_k\right]=R\left[t_k\right]{\left(R\left[t_k\right]\right)}^T$$

(40)

The square root of the state prediction covariance can be defined as follows by introducing the QR decomposition

$$R^{-}\left[t_k\right]={\left(qr\left\{{\left[FR\left[t_k\right]\sqrt{Q\left[t_k\right]}\right]}^T\right\}\right)}^T$$

(41)

In the SRUKF algorithm for 5G NLOS mitigation observations, process noise mainly consists of two parts: deterministic noise and stochastic noise. Deterministic noise can be addressed through prior knowledge, while stochastic noise caused by modeling errors should be accurately compensated. Due to the inaccuracies in modeling, the expected mean and covariance of the innovations may differ from the actual values. We propose a stability coefficient to address this issue.

$$R\left[t_k\right]=\sqrt{\phi \left[t_k\right]}{R}^{-}\left[t_k\right]$$

(42)

The estimated innovation covariance $P_{\gamma ,\gamma }\left[t_k\right]$ should be larger than or equal to the actual one $\mathrm{E}\left\{\gamma \left[t_k\right]{\gamma }^T\left[t_k\right]\right\}$ in order to maintain consistent estimation.

$$P_{\gamma ,\gamma }\left[t_k\right]\ge \mathrm{E}\left\{\gamma \left[t_k\right]{\gamma }^T\left[t_k\right]\right\}$$

(43)

$$\gamma \left[t_k\right]=y\left[t_k\right]-\widehat{y}\left[t_k\right]$$

(44)

The innovation covariance has the following expression:

$$P_{\gamma ,\gamma }\left[t_k\right]=H\left[t_k\right]P\left[t_k\right]H^T\left[t_k\right]+Q\left[t_k\right]$$

(45)

Where $P\left[t_k\right]$ denotes the covariance of the predicted state and the crosscovariance between the measurements and the predicted state is $P_{x,y}\left[t_k\right]$

$$P_{x,z}\left[t_k\right]=P\left[t_k\right]H^T\left[t_k\right]$$

(46)

Combining (45) and (46) leads to (47)

$$P_{\gamma ,\gamma }\left[t_k\right]=P_{x,y}^T\left[t_k\right]P^{-1}\left[t_k\right]P_{x,y}\left[t_k\right]+Q\left[t_k\right]$$

(47)

Incorporating across (43) to (47) yields (48), the fading factor introduced into the predicted state covariance should meets the requirements

$$\phi \left[t_k\right]\ge \frac{E\left\{y\left[t_k\right]y^T\left[t_k\right]-Q\left[t_k\right]\right.}{{\left(P_{x,y}^{-}\left[t_k\right]\right)}^T{\left(P^{-}\left[t_k\right]\right)}^{-1}{P}_{x,y}^{-}\left[t_k\right]}$$

(48)

Where $P^{-}\left[t_k\right]$ represents the covariance of the predicted state $P\left[t_k\right]=$$\phi \left[t_k\right]P^{-}\left[t_k\right],P_{x,y}^{-}\left[t_k\right]$ denotes the cross-covariance of the state and measurement without the fading factor, given by $P_{x,y}\left[t_k\right]=\phi \left[t_k\right]P_{x,y}^{-}\left[t_k\right]$.

As a result, it can be confirmed that the lower bound of the stabilized coefficient $\phi \left[t_k\right]$, as obtained in Equation (48), can be utilized to mitigate modeling errors while maintaining the estimation consistency of the innovation. The square root of the predicted measurement covariance can be derived according to Equation (41).

$$R_{y,y}\left[t_k\right]={\left(qr\left\{{\left[\begin{array}{cc}{y}^{*}\left[t_k\right]\hfill & \sqrt{U\left[t_k\right]}\hfill \end{array}\right]}^T\right\}\right)}^T$$

(49)

$$y^{*}\left[t_k\right]=\left[{\left(y\left[t_k\right]\right)}_0-\widehat{y}\left[t_k\right]\cdots {\left(y\left[t_k\right]\right)}_i-\widehat{y}\left[t_k\right]\cdots {\left(y\left[t_k\right]\right)}_{2L}-\widehat{y}\left[t_k\right]\right]\times diag\left(\sqrt{W^{\left(c\right)}}\right)$$

(50)

Here, $\mathrm{L}$ represents the dimension of the receiver. ${\left(y\left[t_k\right]\right)}_i$, where i ranges from 0 to $2\mathrm{L}$, signifies the propagated sigma point, and $y\left[t_k\right]$ denotes the predicted measurement results. $U\left[t_k\right]$ represents the Gaussian measurement noise matrix at time $t_k$. The term $\sqrt{W^{\left(c\right)}}$ represents the weight of the sigma points. The state covariance derived from QR decomposition can be calculated by:

$$R\left[t_k\right]={\left(qr\left\{{\left[{\chi }^{*}\left[t_k\right]-y^{*}\left[t_k\right]K\left[t_k\right]\sqrt{U\left[t_k\right]}\right]}^T\right\}\right)}^T$$

(51)

$${\chi }^{*}\left[t_k\right]=\left[{\left(\chi \left[t_k\right]\right)}_0-\widehat{s}\left[t_k\right]\cdots {\left(\chi \left[t_k\right]\right)}_i-\widehat{s}\left[t_k\right]\cdots {\left(\chi \left[t_k\right]\right)}_{2n_x}-\widehat{s}\left[t_k\right]\right]\times diag\left(\sqrt{W^{\left(c\right)}}\right)$$

(52)

Where ${\left(\chi \left[t_k\right]\right)}_i$, and i ranges from 0 to $2\mathrm{L}$, denotes the propagated sigma points, $K\left[t_k\right]$ represents filtering gain.

4. Simulation and Experiment Results Analysis

In this section, we evaluate the performance of the proposed NM-SLE algorithm through simulation and field experiments. The effectiveness of the proposed method for NLOS observation identification and suppression is validated in real environments. The algorithm proposed can adequately exploit the features of NLOS signals and effectively suppress the residual NLOS measurement errors in the position estimation process. In both simulation and field experiments, we have compared the performance of the proposed localization algorithm with four existing advanced NLOS error elimination algorithms: SVM-REM[16], Vector Tracking Loop (VTL)[17], Virtual Base station(VBS) [19] and Improved Robust Unscented Kalman filter (IRUKF)[26]. The data set used in the simulation is composed of 12000 LOS samples and 5200 NLOS samples.

4.1. Simulation Results Analysis

In Figure 3, the correlation of candidate features selected for NLOS observations in this paper is presented using a correlation coefficient matrix plot. The correlation matrix in the figure is symmetric about the diagonal, and each color within the square corresponds to the same meaning and value. From the analysis of the figure, it can be observed that there is a significant positive correlation between MPV and TE, MPV and SNR, and MPV and KUR, with MPV showing particularly strong correlation with TE. On the other hand, there is a noticeable negative correlation between RT and SNR, RT and MPV, and RDS and SNR. The correlation among other features is relatively weak.

In Figure 4, the overall performance of different methods are compared in terms of MSE. From the graph, it can be observed that SVM-REM exhibits a significant performance drop as the number of features increases, while the VBS method performs the best. The performance of the NM-SLE algorithm proposed in this paper is close to that of the VBS algorithm, with both methods showing lower errors compared to VTL and IRUKF algorithms. The average positioning errors for the proposed method, SVM-REM, VTL, VBS, and IRUKF algorithms are 1.56m, 8.16m, 1.91m, 1.52m, and 1.85m respectively.

Figure 5 presents the cumulative distribution function (CDF) of positioning errors for the NM-SLE method proposed in this paper, as well as for SVM-REM, VTL, VBS, and IRUKF methods. It can be observed that the performance of the NM-SLE method and VBS is significantly better than the other three methods. At a 60% probability level, both NM-SLE and VBS achieve positioning accuracy better than 2 meters. The 90% positioning accuracy for NM-SLE, SVM-REM, VTL, VBS, and IRUKF methods are 3.31m, 4.59m, 4.57m, 3.42m, and 4.87m respectively.

The average execution time of different NLOS suppression methods is presented in Figure 6. The data shown in the graph were obtained by running the methods on a laptop. Although the absolute values may not be directly applicable to resource constrained embedded systems, the relative performance differences between different NLOS suppression methods still reflect their varying computational complexities.

Although this section has demonstrated the effectiveness of the proposed NM-SLE method through simulations, there are still some limitations that need further investigation. Firstly, due to hardware limitations, the simulations only considered data obtained in static scenarios. Further validation is required in dynamic scenarios. Additionally, different shapes and materials may have significant effects.

4.2. Experiment Results Analysis

In this section, we built an indoor NLOS 5G positioning test platform based on 5G positioning signal transmitters, clock synchronization modules, atomic clocks, and 5G positioning terminals. The transmitters shared the same atomic clock and achieved time synchronization through fiber optic synchronization modules. The positioning terminals measured TDOA values and transmitted them back to a laptop for position calculation. In Figure 7, the left side shows the physical structure of the 5G positioning signal transmitter used in the experiment, which consists of a power supply module, power amplification module, and positioning signal generation module. The right side shows the structure of the 5G positioning terminal. The terminal performs baseband signal processing and data demodulation in the FPGA, algorithm control and position calculation in the ARM processor, and information transmission via the Bluetooth module.

Figure 8 shows the experimental test scene. A total of six 5G positioning signal transmitters are deployed. There are 4 pillars in the center of the test area as wireless signal obstructions. The blue circle in the figure indicates the location of the positioning signal transmitter, and the yellow square indicates the location of the test point where the positioning terminal is placed.

In Figure 9 the CDF of position error is displayed. The performance of each algorithm in the actual environment is slightly lower compared to the simulation results. The NM-SLE algorithm proposed in this paper outperforms the other four algorithms significantly. The positioning error is below 4 meters for over 90% of the cases. Additionally, the theoretical average values of the positioning error calculated from the derived CDF are presented in Table 2. Compared to the other four algorithms, the proposed NM-SLE algorithm exhibits the best performance in field experiment.

5. Discussion

In this paper, we propose a method aimed at addressing the challenge of mitigating NLOS errors in 5G indoor localization. Our approach involves several key steps for 5G terminals: firstly, we perform classification and identification of NLOS observations based on decoupling wireless signal features, where indoor propagation channels are categorized into multiple classes to assess the severity of NLOS effects based on channel identification results. Secondly, we establish a 5G NLOS measurement identification and classification criteria for filtering ranging results. Finally, we introduce a position estimation method based on SRUKF to ensure the numerical stability of filter position estimation output results during the fusion of observations. The theoretical simulations and field experiment demonstrated the superior performance compared to four advanced methods: SVM-REM, VTL, VBS and IRUKF. The proposed method can provide crucial support for high precision localization in commercial 5G networks.