Tensor-Based Subspace Tracking for Time-Delay Estimation in GNSS Multi-Antenna Receivers

Although Global Navigation Satellite Systems (GNSS) receivers currently achieve high accuracy when processing their geographic location under line of sight (LOS), multipath interference and noise degrades the accuracy considerably. In order to mitigate multipath interference, receivers based on multiple antennas became the focus of research and technological development. In this context, tensor-based approaches based on Parallel Factor Analysis (PARAFAC) models have been proposed in the literature, providing optimum performance. State-of-the-art techniques for antenna array based GNSS receivers compute singular value decomposition (SVD) for each new sample, implying into a high computational complexity, being, therefore, prohibitive for real-time applications. Therefore, in order to reduce the computational complexity of the parameter estimates, subspace tracking algorithms are essential. In this work, we propose a tensor-based subspace tracking framework to reduce the overall computational complexity of the highly accurate tensor-based time-delay estimation process.


Introduction
Global Navigation Satellite Systems (GNSS) provide geospatial positioning at any point in the gobe through the use of artificial satellites. These systems allow receptors on the surface of the Earth to determine geographic location by means of time of the transmitted signal time-delay, which is the interval that the signal takes to travel from the satellite to the receiver [1]. There are a wide number of applications for GNSS systems, e.g., civil aviation, navigation, geographic mapping, agricultural research, road flow, and defense. Additional recent applications of GNSS systems include: autonomous vehicles that require high safety and precision standards [2], real-time fishing vessel location for sustainable management of fisheries [3] and precision agriculture [4]. In the latter, the expensive agricultural machines are coordinated to work 24 hours a day and seven days a week and the costs with fertilization and the impact on the environment are drastically reduced due to the precise application of fertilizers.
The accurate positioning provided by GNSS receivers depend on its line of sight (LOS) with the satellites. In practice, non-LOS (NLOS) components are also present due to reflections on buildings, trees, and poles. These multipath components can drastically degrade positioning accuracy. In order

Notation
The notation used throughout this paper is organized as follows: scalars are denoted by lower-case letters: {a, b, c, . . . }, vectors are written as bold lower case letters: {a, b, c, . . . }, matrices by uppercase bold letters {A, B, C, . . . }, and tensors are represented as calligraphic letters: {A, B, C, . . . }. The superscripts: T , * , H , −1 and + , denote the transpose, conjugate, conjugate transpose (Hermitian), inverse, and pseudo-inverse of a matrix, respectively. For a matrix A ∈ C M×N , its m-th row is denoted by A (m , :) and its n-th column is denoted by A (: , n) . The identity matrix is referred by the letter I. The Kronecker product is represented by ⊗ and the Khatri-Rao product (also known as column-wise Kronecker product) by . The Khatri-Rao product between matrices A ∈ C I×L and B ∈ C K×L is a matrix of dimensions IK × L. This product is defined as: The n-th mode unfolding of the tensor A is denoted by [A] (n) [27]. The n-mode product between a tensor A and a matrix B is denoted by A × (n) B.

Data Model
This section is divided in two parts. In Section 2.1 the tensor-based pre-correlation data model is explained. In Section 2.2 the post-correlation signal model constituted by the multiplication of the received signal by a compressed correlator bank is detailed. The data model used in this work is based on [10]. It is important to mention that in the scenario under consideration the model order L d should be known. Note that the estimation of the model order can be obtained using the approaches [28][29][30].

Pre-Correlation Signal Model
In the pre-correlation stage the L received signals for the Uniform Linear Array (ULA)-based GNSS receiver, composed of a LOS component and (L − 1) delayed NLOS replicas, are temporally grouped into K epochs, for k = [1 , 2 , . . . , K]. In each epoch, N samples are collected by the M elements of the array. Therefore, according to [10] the received signal can be modeled as: where the matrixΓ T = γ 1 · · · γ ∈ C K×L d contains the complex amplitudes related to each signal component of the d-th satellite along K epochs with γ = γ [1] · · · γ[K] T ∈ C K being the vector with the complex amplitudes related to each epoch. The matrixÃ = a(φ 1 ) · · · a(φ ) ∈ C M×L d denotes the steering matrix of the k-th observation period of the d-th satellite where each a(φ ) ∈ C M is referred as the steering vector with azimuth angle (φ ). The matrixC T = ∈ R N×L d contains the sampled and delayed PRN sequences of the d-th satellite. The term I 3,L d ∈ C L d ×L d ×L d is the third-order identity tensor and N ∈ C K×N×M refers to the AWGN characteristics of the channel which can be expressed as a tensor.

Post-Correlation Signal Model
As GNSS receivers perform a series of cross-correlations in order to align the incoming C/A code with a local generated replica, a bank of compressed correlators Q (d) [10,31], is employed to multiply the received signal with all possible shifted replicas. This operation computes the cross-correlation vector necessary to perform the time-delay estimation.
The application of the compressed correlator bank on the pre-processed signal in (1) leads to the following post-processed signal model:

Proposed Tensor-Based Subspace Tracking Framework
In this section, the proposed subspace tracking techniques for time-delay estimation are presented. First, the PAST-HOSVD/TDE is derived. Then, the PAST-DoA/KRF is formulated. Finally, computational complexity issues are discussed.
As depicted in Figure 1, the adaptive tensor setting [26] which is also referred as online or incremental setting, is defined as a sequence of R observed tensors represented by Y  In accordance with (3) the expression for each tensor of the sequence Y (d) (r) can be denoted as: In order to apply tensor-based subspace tracking methods, for the sake of simplicity only one epoch is considered in (4), implying into a matrix representation of the data associated to the reverse cyclical third-mode unfolding of Y (d) (r) expressed as:

PAST-HOSVD/TDE
According to [10], the HOSVD-based time-delay estimation (HOSVD/TDE) exploits the multidimensional structure of the received signal by decomposing it into principal singular vectors for each dimension [32]. Since the HOSVD/TDE is not an adaptive approach, a full HOSVD is computed for each new collected epoch, implying into a high computational complexity. In contrast to the state-of-the-art HOSVD-TDE, the proposed Projection Approximation Subspace Tracking (PAST)-HOSVD/TDE updates the multidimensional eigenspace by tracking a projection matrices of the different dimensions, allowing a considerable reduction of the computational complexity with similar accuracy.
The approach starts by separating the highly correlated signal components on each acquired snapshot. Therefore, the pre-processing techniques Forward-Backward Averaging (FBA) [7] and Expanded Spatial Smoothing (ESPS) [8] are applied on the reshaped received signal snapshot matrix [Y (r) ] (3) ∈ C M×Q . These techniques exploit the symmetry of the antenna array in order to further improve the parameter estimation.
In order to apply FBA, two matrices Π M ∈ C M×M and Π Q ∈ C Q×Q are defined. The technique is applied according to the following equation: The ESPS method is employed by means of selection matrices that generate L s subarrays with M s = M − L s + 1 antenna elements. Each selection matrix is defined as: where s = {1, · · · , L s }. Then SPS is applied to the FBA pre-possessed signal shown in Equation (6) according to the following method: By applying the FBA and ESPS, the original antenna array dimension is split into three dimensions, the first one related to the reduced antenna arrays (M s ), the second related to the subarrays (L s ) and the third related to the FBA doubling the array aperture, implying into a dimension of size 2 [33,34]. Note that, bringing back the epoch dimension, W FBA+SPS (r) in (8) can be reshaped into a fourth order tensor Z FBA+ESPS (r) ∈ C 2×Q×M s ×L s which is referred to as the data tensor observed at each new snapshot [35]. Notice that Z FBA+ESPS (r) is a pre-processed and rearranged version of [Y (r) ] (3) ∈ C M×Q in (5). For futher information about ESPS and FBA applied to GNSS, we refer to [36].
Subsequently, the multidimensional filtering proposed in [10] is performed on Z FBA+ESPS (r), which HOSVD is given by: where G ∈ C 2×Q×M s ×L s is defined as the core tensor and U (1) ∈ C 2×2 , U (2) ∈ C Q×Q , U (3) ∈ C M s ×M s and U (4) ∈ C L s ×L s are unitary matrices containing the respective singular vectors u (n) of each n-mode unfolding of each mode of Z FBA+ESPS (r). Since the signal model states that the LOS component has the greatest power when compared with the NLOS replicas, the n-mode singular vectors u (n) are ordered in the unitary matrices U (n) in a decreasing order of the magnitude of its corresponding singular values. Hence, the dominant singular vectors of the first, third and fourth modes of Z FBA+ESPS (r) in Equation (9) will be mostly correlated to the LOS component. The eigenfiltering process proposed in [10] is performed according to: where the term ΣV H ∈ C Q×Q , which multiplies the filter, compensates the correlation process with the bank shown in Equation (3). The output represents the multi-dimensional filtered cross-correlation values at each tap of the correlator bank. Later, the amount of points of the output of the correlator bank is increased via a cubic spline interpolation and one dimension peak search is performed in order to locate the time-delay of the greatest power component.
In order to obtain an estimation of each n-mode singular vector u (n) at time instant (t + 1), PAST-HOSVD/TDE employs the tensor-based version of PAST proposed by [21]. The difference lies in the fact that the subspace tracking algorithm will update the n-mode singular vectors by tracking a projection matrix for the first, third and fourth modes of Z FBA+ESPS (r + 1) ∈ C 2K×Q×M s ×L s . The pseudo-code of the proposed technique can be consulted in Algorithm 1. The matricesÛ (1) S (0) = I L s ×L refer to the signal subspace basis related to the respective first, third and fourth modes of Z FBA+ESPS (r) in (9) andĈ (1) yy (0) =Ĉ (3) yy (0) =Ĉ (4) yy (0) = I L is a d × d correlation matrix of the projection matrix Y(t) which is recursively updated in order to track the signal subspace associated to the corresponding unfolding of Z FBA+ESPS (r) at a given time instant (t + 1). After obtaining an estimation for the signal subspace basis, the dominant eigenvectors are selected to perform the eigenfiltering process seen in (10). Note that the tensor-based version of PAST can run on the unfoldings of Z FBA+ESPS (r) in parallel, which provides an efficient implementation.

PAST-DoA/KRF
Since the post-correlation received signal model in Equation (1) follows a PARAFAC model, DoA/KRF proposed in [11] first estimates the factor matrices and then uses one factor matrix to estimate the time-delay. The approach starts by reconstructing the array steering factor matrix by means of a DoA estimation method referred as the Estimation of Signal Parameters via Rotational Invariance Techniques (ESPRIT) [37]. The remaining factor matrices are then obtained by means of the Khatri-Rao factorization (KRF). Later, the time-delay estimation is performed via the signal's greatest power approach, as discussed on the previous section. The pseudo-code of DoA/KRF is provided in Algorithm A1.
When considering an adaptive tensor setting for real-time applications, DoA/KRF [11] updates the signal subspace estimation by performing a truncated SVD on each observed epoch or on a K collection period of epochs. In this context, to provide a reduction on the overall computational complexity, the proposed PAST-DoA/KRF method links the Tensor-based PAST [35] algorithm to the DoA/KRF time-delay estimation technique by providing a subspace basis to ESPRIT at each sampling instant. The pseudo-code of the proposed technique can be consulted in Algorithm 2. Notice that the tensor-based PAST technique consists of a rank-1 update of the estimated signal subspace matrix. The method does not provide the exact basis of eigenvectors related to the signal subspace of the covariance matrix of the observed data, but a mere basis to span it. Hence, since any basis belonging to span(U s ) is a solution for the shift invariance equation of ESPRIT [37], the array steering factor matrix can be rebuild by extracting information about given signal subspace estimation related to the antenna array dimension. We refer to [38] for more detailed information about ESPRIT applied to the estimated factor matrices in order to estimate the DoA.

Computational Complexity
In this assessment, the operation counts on the study of the computational complexity of the algorithms are expressed in terms of Multiply Accumulate (MAC) operations, referred as FLoating point OPeration (FLOP) counts [39]. Similarly as considered in [10], the objective of the computational complexity evaluation is to provide an assessment about the relative complexity difference between each algorithm. Since the proposed techniques differ from the state-of-the-art only when it concerns the update of the signal subspace estimation, operations such as the unfolding and the inverse-unfolding are not taken into account. In addition, the stage of finishing the correlation by multiplying the solution vector by ΣV H ∈ C Q×Q will also be ignored once it is an operation performed by all algorithms. For additional information concerning the absolute computational complexity of the state-of-the-art algorithms, we refer to [12].
It can be noticed from Table 1 that if the receiver decides to increase the number K of collected epochs, the computational cost of the proposed methods remains unaltered while the complexity of the state-of-the-art techniques increase.
% -Tracking the signal subspace related to the dimension M s .
% -Tracking the signal subspace related to the dimension L s .

Simulations
Similarly to [10,11] we consider a scenario with a left center-hermitian uniform linear array receiver with M = 8 elements and half-wavelength ∆ = λ/2 spacing. Note that the GALANT hardware proposed by the German Aerospace Center (DLR) in [40] allows up to 16 antennas. Typical range goes from 4 to 6 antennas. The GNSS signal is a GPS second generation course acquisition PRN code from d = 1 satellite at a carrier frequency f c = 1575.42 MHz, bandwidth B = 1.023 MHz and chip duration T c = 1/B = 977.52 ns with N = 2046 samples collected every k-th observation period during K = 30 epochs. Each epoch has duration ∆t = 1 ms. The number of impinging signals on the receiver is L = 2, one LOS component with time-delay τ (LOS ) and one NLOS multipath replica with time-delay τ (NLOS) , such that τ (NLOS) = τ (LOS) + ∆τ, where ∆τ is the delay difference between each component. Their azimuth angle difference is ∆φ = 30 • . For the pre-processing techniques SPS/ESPS the array is divided into L s = 5 subarrays with M s = 4 elements each. Signal phases are independent and identically distributed ∼U[0, 2π[. The number of correlators in the bank is Q = 11 equally spaced between −T c and T c . The carrier-to-noise ratio is C/N 0 = 48 dB-Hz, resulting in a pre-correlation SNR pre = C/N 0 − 10 log 10 (2B) ≈ −15.11 dB. Post-correlation SNR post = SNR pre + G ≈ 15 dB and signal to multipath ratio (SMR) of 5 dB. For the tracking scenario, we consider a total of N = 50 received tensors. The step variation of the angle of arrival between each acquired tensor is ∆φ n = 0.25 • . For the Tensor-based PAST subspace tracker, the forgetting factor β PAST is set to 0.95.
The results are obtained by performing M c = 1000 Monte Carlo (MC) simulations to plot the root mean squared error (RMSE) expressed in meters, which is a measurement of the estimated time-delay multiplied by the speed of light, c = 299792458 m/s. The expression is written as: The RMSE (τ LOS ) [m] for the techniques DoA/KRF, PAST-DoA/KRF, HOSVD/TDE and PAST-HOSVD/TDE simulated in the tracking scenario with the parameters described before is presented in Figure 2    In order to evaluate the subspace estimates provided by the proposed techniques, Figure 4 displays the Largest Principle Angle (LPA) between the true subspace and the estimated signal subspace along the number of observed tensors. The LPA between the column spaces of two matrices U 1 , U 2 ∈ C M×d is computed based on [39] according to: where orth{·} is the orthonormal basis operator and σ min {·} is the smallest singular value operator.
Notice that the LPA of PAST-DoA/KRF is smaller than the LPA of the estimated basis of PAST-HOSVD/TDE related to the modes of the observed tensor. This difference reflects on the performance of the RMSE∼[m] as observed in Figure 2.

Conclusions
In this work a tensor-based subspace tracking framework for time-delay estimation in GNSS receivers is proposed. The proposed PAST-DoA/KRF and the state-of-the-art DoA/KRF present similar accuracy, although the computational complexity of the proposed approach is considerably smaller than the state-of-the-art counterpart. Nonetheless, when comparing the performances of PAST-HOSVD/TDE and HOSVD/TDE it is observed that a small positioning error is produced. This difference relies on the fact that in the HOSVD-based approach a larger number of subspaces are tracked. On the other hand, on the computational complexity of the simulation, the reduction of the cost is more expressive.
The proposed methods reduce the cost, allowing best performance for applications that require real-time processing. Note that in this work only the additive noise and the multipath components are considered. The proposed techniques need to be adapted when taking into account propagation effects such as refraction in the ionosphere or troposphere. These and other related effects are considered as future work.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this paper:

Appendix A
This appendix contains the pseudo-code for DoA/KRF [11] which can be consulted in Algorithm A1.
1. ESPRIT [37] section: Estimate the angle of arrival of the L received signals.
Initialization: 1.1 -Define the selection matrices J 1 and J 2 .
Obtain an estimation for the signal subspace. U s = EVD (R xx ) Solve the shift invariance equation to obtain an estimation for Ψ. Ψ = (J 1Ûs ) + J 2Ûs Obtain an estimation for the spatial frequenciesμ . {μ } = diag EIG Ψ Find the corresponding angles of arrival {θ }.