Tensor-Based Framework With Model Order Selection and High Accuracy Factor Decomposition for Time-Delay Estimation in Dynamic Multipath Scenarios

Global Navigation Satellite Systems (GNSS) are crucial for applications that demand very accurate positioning. Tensor-based time-delay estimation methods, such as CPD-GEVD, DoA/KRF, and SECSI, combined with the GPS3 L1C signal, are capable of, significantly, mitigating the positioning degradation caused by multipath components. However, even though these schemes require an estimated model order, they assume that the number of multipath components is constant. In GNSS applications, the number of multipath components is time-varying in dynamic scenarios. Thus, in this paper, we propose a tensor-based framework with model order selection and high accuracy factor decomposition for time-delay estimation in dynamic multipath scenarios. Our proposed approach exploits the estimates of the model order for each slice by grouping the data tensor slices into sub-tensors to provide high accuracy factor decomposition. We further enhance the proposed approach by incorporating the tensor-based Multiple Denoising (MuDe).


I. INTRODUCTION
As Global Navigation Satellite Systems (GNSS) become more ubiquitous and this technology proved to be essential The associate editor coordinating the review of this manuscript and approving it for publication was Wei Wang .
for applications such as civilian aviation, autonomous driving, defense, and timing and synchronization of critical networks. GNSS receivers require line of sight (LOS) signals from at least four satellites to estimate their position on the Earth's surface. Additionally, besides the LOS component, non-LOS (NLOS) multipath components are received VOLUME 8, 2020 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ due to the reflections on trees, poles, lamps, and buildings. Therefore, the superposition of the LOS and NLOS multipath components degrades the time-delay estimation (TDE) and, consequently, the positioning estimation. State-of-the-art GNSS receivers equipped with a single antenna are in general remarkably sensitive to the effect of multipath components [1]- [3]. Thus, multi-antenna GNSS receivers became the focus of research on resilient positioning withstanding not only multipath but also interference and spoofing. In addition to beamforming approaches [4], [5], multi-dimensional parameter estimation approaches [6], [7], and other approaches as those proposed in [8] and [9], tensor-based decomposition methods showed significant improvements over matrix-based decomposition methods. A tensor-based decomposition features uniqueness, improves the identifiability of the parameters, and the tensor structure permits efficient denoising of the received signal. Therefore, tensor-based multipath mitigation methods, combined with antenna arrays, have been proposed as an alternative to single antenna and matrix-based techniques.
In [10], the authors propose a so-called Tensor-based Eigenfilter using the Higher-Order Singular Vector Decomposition (HOSVD) combined with Forward-Backward Averaging (FBA) [11], Spatial Smoothing (SPS) [12], [13], and a bank of correlators to mitigate multipath and to improve time-delay estimation of the LOS signal. In [14] a three step approach based on direction of arrival (DoA) estimation, the Khatri-Rao factorization (KRF), and a bank of correlators was proposed. This DoA/KRF method [14] outperforms traditional matrix-based decomposition methods. However, the performance of this approach depends on estimates of the DoAs of the received signals to construct the respective loading matrices of the signal tensor decomposition. Furthermore, [14] assumes a static environment, thus it cannot be applied to scenarios where the model order changes overtime. In [15], the authors proposed the Canonical Polyadic Decomposition by a Generalized Eigenvalue Decomposition (CPD-GEVD). Although the CPD-GEVD showed robustness in the presence of multipath components and array imperfections, it assumes a constant model order between data epochs. In [16], the tensor-based methods proposed in [10], [14], [15] were extended to third-generation GPS (GPS3) signals, since GPS3 is more robust against multipath components in comparison to the signals of the second-generation GPS (GPS2) due to its Time Multiplexed Binary Offset Carrier (TMBOC) modulation [17]- [19]. In [20], the authors showed that the SEmi-algebraic framework for computing an approximate CP decomposition via Simultaneous Matrix Diagonalizations (SECSI) [21] can be applied to GPS2 and GPS3 signals and antenna array-based receivers using GPS3 and SECSI outperform antenna array receivers using GPS2 and SECSI. This underlines the superior multipath performance of the new GPS3 signals. Additionally, [20] showed that the SECSI method is more robust in the presence of strongly correlated signals than the CPD-GEVD method while on the other hand being computationally more expensive.
The Tensor-based Eigenfilter [10] does not require a previous model order selection (MOS), but achieves lower accuracy compared to the methods proposed in [15], [20]. However, the CPD-based techniques, introduced in [15], [20], require an estimate of the tensor model order to enable tensor factorization. Thus, even though [15], [20] have improved the performance of TDE in presence of highly correlated multipath, these methods assume a known and constant model order. Recently, in [22], the authors proposed the (L r , L r , 1)-GEVD approach to perform multi-linear rank-(L r , L r , 1) decomposition by clustering NLOS components and obtaining the LOS component. Meanwhile, in [23], the authors proposed a tensor-based subspace tracking framework to keep track and update the tensor signal subspace. However, similarly to previous tensor-based methods, [22] and [23] assume that the model order is constant between data epochs. Hence, MOS schemes need to be included and assessed for the application of tensor-based methods for GNSS in order to achieve highly accurate and robust TDE of the LOS signal in case correlated multipath signals are received.
The literature on MOS for matrix data models is quite extensive. The following approaches can be considered the state-of-the art of matrix-based MOS: 1-D Akaike's Information Criterion (AIC) [24], 1-D Minimum Description Length (MDL) [24], EFT [25], Modified EFT (M-EFT) [24], RADOI [26], and the subspace-based ESTER [27]. Also, tensor-based MOS was proposed as well, such as R-D AIC [28], R-D MDL [28], and R-D EFT [24]. Even though MOS methods are crucial to mitigate multipath components, the literature hardly investigates MOS methods for GNSS applications. For example, in [29], the authors apply MDL with antenna arrays to estimate GNSS NLOS components, thus showing that the model order estimation depends on the carrier-to-noise density ratio (C/N 0 ) of the received signals. Furthermore, in [30], a multipath detection algorithm is developed based on a one-way analysis of the variance (ANOVA). However, in [30], the authors assume the DOA of the LOS signal to be known.
In this paper, we propose a tensor-based framework with model order selection and high accuracy factor decomposition for TDE in dynamic multipath scenarios for GNSS. This new approach includes the estimation of the model order for each slice of the data tensor and subsequent grouping of the data tensor into sub-tensors. To successively apply SPS, denoising, and reconstruction of the noisy data, the proposed high accuracy tensor-based decomposition, uses the respective sub-tensors combined to the tensor-based MUltiple DEnoising (MuDe) [31]. Therefore, the proposed framework is composed of four steps, namely, estimation of the number of multipath components in slow and fast dynamic multipath scenarios, further mitigation of the multipath effect by employing the tensor-based MuDe, separation of the sources by using our proposed high accuracy tensor-based decomposition exploiting the different model orders of the slices of the data tensor, and estimation of the time-delay of the signal components. This paper has six sections, including this introduction. Section II introduces the notation utilized in the paper. In Section III, the data model assuming multipath components and the dynamic model order is presented. Section IV describes the proposed tensor-based framework with model order selection and multiple denoising for dynamic multipath components. Section V presents and discusses simulation results for the proposed framework. Section VI draws the conclusions.

II. NOTATION
In this section, we introduce the mathematical notation used in this paper. Scalars are represented by italic letters (a, b), vectors by lowercase bold letters (a, b), matrices by uppercase bold letters (A, B), and tensors by uppercase bold calligraphic letters (A, B). The superscripts T , * , H , −1 , and + denote the transpose, conjugate, conjugate transpose (Hermitian), the inverse of a matrix, and pseudo-inverse of a matrix, respectively.
Additionally, for a vector a ∈ C N , the nth element is denoted as a n . Moreover, for a matrix A ∈ C M ×N , the element in the mth row and nth column is denoted by a m,n , its mth row is denoted by A m,· , and its nth column is denoted by A ·,n . Furthermore, the vec{·} operator reshapes a matrix into a vector by stacking all elements in a column vector. The operators and • denote the Khatri-Rao product and the outer product, respectively.
The operator cond{·} computes the condition number of a matrix. The smaller is the condition number, the more is stable the inverse, of the said matrix [32]. The Frobenius norm of a matrix A is denoted by ||A|| F whereas ||A ·,n || 2 is the vector norm. For a matrix A ∈ C M ×N with M < N , the diag{·} operator extracts the diagonal of a matrix.
The n-mode unfolding of tensor A is denoted as [A] (n) , which is the matrix form of A obtained by varying the nth index along the rows and stacking all other indices along the columns of [A] (n) . The n-mode product between tensor A and a matrix B is represented as A × n B. The N th order I N ,L ∈ R L×...L is defined as the N -way identity tensor of size L, whose elements are equal to one if the N indices are equal and zero, otherwise.

III. DATA MODEL
This section firstly introduces the scenario considered in this work in Subsection III-A. Subsection III-B describes how the signal tensor is constructed. Finally, in Subsection III-C, the post-correlation data model is defined for a GPS L1C pilot signal [17].

A. SCENARIO
We consider a GNSS receiver equipped with an antenna array with M elements. We assume that for the received signals of d = 1, . . . , D satellites, the LOS signal of the dth satellite is superimposed with L d (k) − 1 NLOS multipath components. The observations are collected during K periods (or epochs) each with N samples, where k = 1, . . . , K and n = 1, . . . , N . Moreover, the total number of received signal components is L (k) = D d=1 L d (k) , where (k) = 1, . . . , L (k) and d (k) = 1, . . . , L d (k) are the th and d th component at the kth epoch. Furthermore, we assume that τ As shown in [15], a tensor model can be used to express the received complex baseband signal at the output of the M antennas of an antenna array as collects the complex amplitudes of each signal component during K epochs with including the complex amplitudes at each epoch. The matrix (4) comprises the sampled L1C pilot code sequence with ] ∈ C N ×1 collecting the periodically repeated pseudo random binary sequences (PRBSs) with time-delay τ for each satellite and the respective multipath components. Then, the matrixÃ collects the array responses a(φ (k) ) ∈ C M ×1 with azimuth angle φ (k) of (k) components at the kth epoch. Moreover, N is a white Gaussian noise tensor. The tensor in (1) is composed of three dimensions, the first dimension of size K is related to each epoch, the second dimension of size N is associated with the collected samples in each epoch, and the third dimension of size M corresponds to the spatial diversity of the receive antenna array.

C. POST-CORRELATION DATA MODEL
To separate the L d (k) LOS signals and NLOS multipath components of the dth satellite, the GNSS receiver uses a bank of correlators with respect to each satellite. Thus, the GNSS receiver applies D banks of correlators on the received signal, obtaining D output signals. We define the dth bank of correlators with Q ''taps'' as  with τ 1 < . . . < τ Q and the qth delayed reference sequence c d [τ q ] corresponding to the qth tap. The thin SVD of Q d with provides the bank of correlators which performs cross-correlation (compression) preserving the noise statistics [33]. Therefore, Q ω suppresses other satellites while preserving the noise statistics. When comparing Q (d) ω to the bank of correlators Q d , we observe that Q d provides a sampled cross-correlation function of the bank of correlators with the received LOS component. Thus, according to [10], we can correlate the received signal tensor ω to separate the dth satellite from all other received satellites and we obtain where holds the cross-correlation values of the LOS and NLOS components for satellite d, and A ∈ C M ×L d (k) comprises the L d (k) array responses of satellite d excluded fromÃ. In the following we assume that A does not vary significantly over K epochs. N ω ∈ C K ×Q×M denotes a white Gaussian noise tensor after correlation. The tensor M is the multiple access interference (crosscorrelation) of the other satellites and their respective multipath components. In general, M is negligible in comparison to other terms, since signals are more or less decorrelated.

IV. PROPOSED FRAMEWORK
As illustrated in Figure 1, the proposed framework can be divided into six steps. In Subsection IV-A we describe block (1). In this step, we receive the post-correlated signal and executes the M-EFT method to estimate the number of multipath components for the dth satellite in dynamic environments. Additionally, in Subsection IV-B, we describe an alternative method for static scenarios, in which we use the RADOI method to estimate the number of multipath components. Next, in Subsection IV-A, we describe block (2) of the framework. In this step, we use the post-correlation tensor Y and the estimated model order L d (k) to generate sub-tensors. In Subsection IV-C, we detail the denoising step in the block (3) by applying the MuDe technique to filter the sub-tensors. Then, block (4), as described in Subsection IV-D, uses the estimated grouped model order and the sub-tensors, filtered by MuDe, to perform the Mode 1 HOSVD SECSI with left-hand matrix (Mode 1 HOSVD SECSI) method to estimate the factor matrices of each sub-tensor. Afterward, as described in Subsection IV-E, block (5) describes how the factor matrices of each sub-tensor are normalized and used to extract the LOS components of all sub-tensors. Finally, as described in Subsection IV-F, we perform TDE of the LOS component of each sub-tensor, as illustrated in the block (6) of the framework given in Figure 1.

A. MODEL ORDER SELECTION FOR DYNAMIC ENVIRONMENTS
To compute the model order applying the Modified Exponential Fitting Test (M-EFT) [24], we use the covariance matrixR[k] obtained from each epoch of the tensor Y (d) as defined in (9). Therefore, we firstly compute the eigenvalue decomposition (EVD) to obtain is a unitary matrix containing the eigenvectors, = diag{λ 1 , . . . , λ M } ∈ C M ×M is a diagonal matrix including the the sorted eigenvalues λ i , such that λ 1 > λ 2 > · · · > λ M , and the covariance matrix the column space of the steering matrix A span the same subspace. Furthermore, the M-EFT can adopt an exponential profile to approximate the Wishart profile of the noise eigenvalues and consequently enabling their prediction. The M-EFT estimates the model order by computing the distance from λ M −P , calculated from the measurements to the predicted eigenvalueλ M −P , where P is a possible number of noise eigenvalues. Furthermore, the M-EFT method computes the threshold coefficients η P , and then estimates the model order.
Since M , Q, and the probability of false alarm P fa , do not vary, η P is computed previously and stored. To compute the M-EFT we refer to [24], [34]. Then, once we have the estimated model orderL d (k) (k) for each epoch, we group the epochs with the same estimated model order, as illustrated in Figure 2(a) and 2(b). Moreover, we create a vector that contains the grouped model order L d (k) (k) (t) with respect to the epochs and for t = 1, . . . , T sub-tensors. Therefore, we createỸ (t) sub-tensors which will be used to perform TDE. For instance, if we estimate three different model ordersL d 1 (k) = 2 andL d 2 (k) = 3 we concatenate the epochs with the same model order to create a new tensor. Afterwards, we can use the sub-tensors to perform TDE.

B. MODEL ORDER SELECTION FOR STATIC ENVIRONMENTS
In order to decompose the tensor Y (d) into factor matrices for TDE, first, the number of multipath components of the dth satellite L d (k) have to be estimated. To perform MOS in static environments we propose to apply RADOI [26] for which the covariance matrixR yy obtained from the thirdmode unfolding of tensor [Y (d) ] (3) as defined in (9) is used to compute an EVD. Thus, we obtain whereR yy ∈ C M ×M is a Hermitian matrix,ũ = [ũ 1ũ2 . . .ũ M ] ∈ C M ×M is an unitary matrix containing the eigenvectors,˜ = diag{λ 1 , . . . ,λ M } ∈ C M ×M is a diagonal matrix collecting the sorted eigenvaluesλ i , such thatλ 1 >λ 2 > · · · >λ M , and the covariance matrix R R qq ∈ C M ×M of the bank of correlators. Moreover, we definẽ U (s) = [ũ 1ũ2 . . .ũ P ] ∈ C M ×P as the truncated matrix composed of P eigenvectors ofŨ corresponding to the P largest eigenvalues of˜ . Therefore, as discussed above for the EFT, in case that P = L d (k) , the dominant eigenvectors U (s) R ∈ C M ×L d (k) and the column space of the steering matrix A span the same subspace.

C. TENSOR-BASED MULTIPLE DENOISING
Since TDE performance is sensitive to signal-to-noise ratio (SNR) and degrades in noisy scenarios, we include a denoising step. Therefore, we propose to use the MuDe approach [31], which is a pre-processing technique to denoise tensor-like data. MuDe combines the principle of Spatial Smoothing (SPS) [12] with successive SVD-based low-rank approximations of the output signals for sub-arrays of varying size in each spatial dimension of the obtained signal tensor and, then, rebuilds the the different sub-arrays into a tensor.

D. ESTIMATION OF FACTOR MATRICES
Originally, the state-of-the-art SECSI method [21] offers a trade-off between performance and reliability. The authors in [21] prove that for a 3-way tensor, one can construct six distinct diagonalization problems by computing the slices of the core tensor. Then, one can use each core tensor to compute the right-hand and left-hand matrices of each slice. Following the solution of the six distinct diagonalization problems, the authors in [21] propose to analyze the estimates and, then, select the estimate with the lowest error. However, to reduce processing time, [20] proposes to utilize only the right-hand matrix obtained from the first-mode slice of the core tensor. Moreover, in [20] it was shown that the SECSI method combined with HOSVD (HOSVD SECSI) introduced in [21] presents a higher complexity than the CPD-GEVD. However, despite the higher complexity of the HOSVD SECSI, this method shows better performance in scenarios with highly correlated signals. Consequently, the HOSVD SECSI method is more reliable in more demanding scenarios. However, based on simulations, we show that, by applying different tensor modes, we can also achieve better performance in dynamic situations.
Since we create a sub-tensor for the tensor epochs with the same model order, in some situations, we obtain a model order greater than the number of epochs in a given sub-tensor. Since the method used in [20] is no longer suitable for a dynamic scenario, we propose to apply a different tensor mode to perform factor matrix estimation. Hence, in several simulations, the left-hand matrix obtained from the thirdmode slice of the core tensor (Mode 1 HOSVD SECSI) S (d) (t) proved to be a suitable alternative. Finally, to compute the Mode 1 HOSVD SECSI we firstly compute the HOSVD lowrank approximation of the sub-tensors Y (d) (t) obtained from grouping the tensor Y (d) epochs according to their model order.
where the superscript (t) indicates the tth sub-tensor, U the truncated singular matrices related to the tth sub-tensor. Moreover, to compute the singular matrices and the core tensor, the SVD is derived for each of the unfolding of the tensor for each dimension. Thus, we can write are the left singular vector matrix, singular value matrix, and right singular vector matrix for the ith dimension, respectively. Consider that the left singular vector matrices can be sequentially computed as follows: Next, let us consider the first-mode unfolding of Y (d) and use the representation in (9) and (14), then, we obtain Hence, the subspace spanned by the columns of U represent the loading matrices of the CPD of the core tensor of the HOSVD given in (14). Thus, the tensor S can be represented as The transform matrices T (22). In theory, the CPD can be directly applied to Y (d) , to extract the factor matrices. However, as demonstrated in [10], by directly computing the factor matrices using Alternating Least Squares (ALS), there are convergence problems, resulting into poor performance in terms of TDE. Therefore, to perform the CPD, it is sufficient to compute the loading matrices T As a consequence of the symmetry of the SECSI problem, we can build 6 Simultaneous Matrix Diagonalization (SMD) problems for a three-way model [21]. However, as shown by [20], without loss of generality, we can only use one mode of the compressed core tensor S (d) (t) to compute the right-hand and left-hand matrices utilized in the SMD step described in [21]. In [20], the authors show that the SECSI third-mode, of a given compressed core tensor, in the GNSS case, yields the best estimation performance in static scenarios. However, after performing several numerical simulations, we now select the left-hand matrix of the first mode of the compressed core tensor S (d) (t) . Therefore, the ith slice of the first-mode of tensor S (d) (t) is selected to compute the lefthand matrix where e T i is a vector with all zeros except in the ith position. Next, we select the slice of the tensor S (d) (t) with the smallest condition number where p is an arbitrary index between one and the total number of slices to be diagonalized and defines the slice of the tensor S t with the smallest condition number where cond{·} computes the condition number of a matrix. The smaller the condition number of a matrix, the more stable is its inversion. Furthermore, we obtain the left-hand matrices S (t),lhs 1,i by multiplying S 1,i by S 1,p on the left-hand side Since p is fixed, we can i in (29) obtaining N −1 equations, since i = p. Our goal is to findT (d) (t) 3 that simultaneously diagonalizes the N − 1 equations. We refer here to the techniques in [35] and [36].
Since in the noiseless case, according to (9), the third-mode unfolding exposes the factor matrix A (t) , we can write Using from (14) andT 3 from the diagonalization step we obtain Afterwards, we multiply (30) by the pseudo-inverse of the estimatedÂ (t) from the left-hand side to obtain Then, the factor matrices (CQ ) T and (t) can be estimated from (32) by applying the Least Squares Khatri-Rao Factorization (LSKRF) [37].

E. LOS SELECTION
Subsequently, to estimating all the parameters of the received signal, we need to separate the LOS and NLOS signal parameters. In this subsection, we describe the fifth element of the framework. This element performs the LOS selection based on the estimated factor matrices. Therefore, as described in [20], to perform the LOS selection the estimated factor matrices (ĈQ (d) ω )) (t) T ,Â (t) , andˆ (t) are normalized to unit norm for the (t) d th component and can be given as Next, with the normalized factor matrices, we construct the tensor G where where } and we obtain Assuming that the received signal component with the largest power corresponds to the LOS signal, we select the respective column of the estimated (ĈQ where q contains the cross-correlation values at each tap of the correlator bank. Then, as shown in [10], [14], [15], [20], [38], [39], the resulting vector q is interpolated using a simple cubic spline interpolation. Thus, by using the resulting interpolated vector, we can derive the cost function F(κ) [20], which is the cross-correlation function with the received LOS signal. Finally, we use this cost function to estimate the timedelay of the LOS signal by solvinĝ and with a sampling frequency f = 2B MHz. Therefore, N = 245520 samples are collected for the L1C pilot code per epoch. The carrier-to-noise ratio is C/N 0 = 48 dB-Hz, resulting in a pre-correlation signal-to-noise ratio SNR pre = C/N 0 − 10 log 10 (2B) ≈ −25.10 dB for GPS3. Given the processing gain G = 10 log 10 (Bt) ≈ 50.9 dB for GPS3. Hence, the postcorrelation signal-to-noise ratio SNR post = SNR pre + G ≈ 25 dB. Moreover, the signal-to-multipath ratio SMR 1 = 5 dB for L d = 2. In case L d = 3 the SMR 1 = 5 dB for the first NLOS signal, and a SMR 2 = 10 dB for the second NLOS signal. Besides the simulation considering exact knowledge of the array response (perfect array), we added errors in the antenna array geometry to distort the real array response with respect to the known response by displacing the antennas in x and y positions according to a normal distribution ∼N (0, σ 2 ). The standard deviation is computed in terms of the probability p = P(e > λ/2), where the error exceeds a half wavelength. We fix the relative time-delay τ at 0.5T c while varying the error probability p from 10 −6 to 10 −1 . Moreover, we consider a probability of false alarm P fa = 10 −3 . Additionally, we performed simulations considering a dynamic scenario with one satellite with PRBS = 17. In this dynamic scenario, we vary the DoA of the LOS component for each epoch and the number of LOS and NLOS components within the tensor Y. We define a DoA difference between epochs of 2 • . Moreover, we define the first 15 epochs have L d = 5 while the last fifteen collected epochs have L d = 6. Since we consider a dynamic scenario, we target a lower probability of false alarm P fa = 10 −6 .
Since [20] showed that simulations potentially have outliers in case the signals are strongly correlated, we, therefore, performed 1000 Monte Carlo (MC) simulations to compare all approaches in terms of the Root Median-Squared Error (RMDSE) of the time-delay of the LOS component. The TDE performance of the proposed tensor-based methods is compared to the TDE performance in the case A and are considered known, which can be considered as a lower bound for TDE performance for the proposed methods.

A. PROBABILITY OF DETECTION CONSIDERING A STATIC SCENARIO
In this section, we present the probability of detection (PoD) for simulations considering an array of antennas for which the array response is known. We compute the PoD for a relative delay τ = 0.5T c . In TABLE 1, we show the PoD for static scenarios with d = 1 satellite with L d = 3 components and d = 3 satellites with L d = 3 components. Note that, in the simplest scenario, d = 1 and L d = 3, only AIC presents a PoD below 99%. However, observe that if we add more satellites to the simulations, the EFT, AIC, and MDL-based methods do not correctly estimate the model order. However, note that both RADOI and ESTER methods have a consistent performance achieving the same PoD for all simulations.

1) PoD CONSIDERING AN ARRAY WITH ERRORS AND A STATIC SCENARIO
In this section we consider an antenna array with errors in its array response model (imperfect array) with d = 1 and L d = 2, d = 1 and L d = 3 impinging signals, and a fixed relative time-delay τ = 0.5T c . In TABLE 2, we show the PoD for the MOS methods for d = 1 satellite with L d = 3 impinging signals and d = 3 satellite with L d = 3 impinging signals.
Note that the ESTER method shows an inferior performance when we add a second NLOS component to the simulation. Observe that in case the array response model of the antenna array includes errors (e.g. the antenna elements of the array experienced displacement errors), the eigenvalue-based methods EFT, M-EFT, MDL, R-D AIC, R-D MDL, R-D EFET, and RADOI, are insensitive to these errors in the array response model. However, since the ESTER method assumes the matrix A has a Vandermonde structure, errors in the array response model cause that the matrix A has a different structure than Vandermonde. Therefore, the ESTER shows to be sensitive to array imperfections. Moreover, we can observe that RADOI is the most accurate MOS method in a static environment.

B. TDE CONSIDERING A STATIC SCENARIO
In this section, we present simulation results of the various methods for the GPS3 L1C signal. We assess the case L d = 3 and d = 3. Additionally, preceding the time-delay estimation, we apply the RADOI method since other matrix and tensor-based methods, do not provide reliable estimates. In Figure 3 we show the TDE error for matrix-based MOS methods in a static scenario. Note that the results obtained by RADOI present similar estimates to the known model order case. Moreover, in the case signals are weakly correlated, e. g., τ > 0.2Tc, the RADOI method provides precise model order estimates. Note that in a static scenario with constant model order, the MOS accuracy seems irrelevant since the TDE performance is similar. The larger power of the LOS component after signal correlation explains the similar TDE. However, the Mode 1 HOSVD SECSI shows higher estimation errors when combined with the AIC since AIC provides an estimated model order of L d > 5. Thus, to keep the TDE error as low as possible, it is essential to use a reliable MOS scheme.

C. PROBABILITY OF DETECTION CONSIDERING A DYNAMIC SCENARIO
In this section, we present the PoD computed considering a dynamic scenario and a perfectly aligned array of antennas with no errors in the array response model. In dynamic scenarios, the number of LOS and NLOS components within the tensor Y (d) are changing. When applying matrix-based MOS methods, we performed the MOS for each epoch and individually computed the PoD. Moreover, in a dynamic scenario, we use the pre-processing methods FBA, SPS, and their combinations to attempt to improve the model order estimation accuracy in case of highly correlated signal components.
In Figure 4 we show the PoD for the M-EFT method for τ = 0.1 T c and K = 30. Note that the M-EFT method  fails to estimate the model order but applying FBA and SPS, M-EFT shows a PoD above 20% when L d = 5 and above 50% when L d = 6.
In Figure 5, we show the M eigenvalues and the M predicted eigenvalues, at τ = 0.1 T c . The eigenvalues are normalized with respect to the strongest eigenvalue (LOS eigenvalue). In Figure 5 one can observe, that FBA and SPS add some gain to the NLOS components. However, since the NLOS components' power is too small, we obtain low NLOS eigenvalues. Note that the 3rd, 4th, and 5th eigenvalues are as weak as the noise eigenvalues, i.e., 6th, 7th, and 8th eigenvalues. Although the eigenvalues are weak, note that, when we combine the FBA and SPS, we obtain a significant gain to allow a more accurate MOS. Therefore, the AIC, ESTER, and RADOI methods are not suitable to perform MOS for GPS3 signals in dynamic environments. MDL, EFT, and M-EFT also show poor performance. However, when applying the pre-processing methods FBA and SPS, we can improve the MOS performance. Mainly, we can achieve better results when combining these pre-processing methods with the MDL, EFT, and M-EFT methods. Since the matrix-based MOS methods use slices of the full tensor, the additional simulations showed that the PoD remains constant as we increase the number of collected epochs. Moreover, we performed simulations with τ > 0.1T c , and we observed that the PoD increases as we increased τ . Therefore, in case the LOS and NLOS components are separated adequately in time, the MOS methods will show better performance, as the different signal components are less correlated in time.

D. TDE CONSIDERING A DYNAMIC SCENARIO
In this section, we present the TDE performance for simulations considering a dynamic scenario and a perfectly aligned array of antennas with no errors in the array response model.
Since the matrix-based AIC, ESTER, and RADOI methods showed poor performance in a dynamic scenario, we only used the matrix-based MDL+FBA+SPS, EFT+FBA+SPS, and EFT+FBA+SPS methods. Additionally, we performed simulations considering the model order to be known for each epoch. Thus, we could divide the main tensor into new sub-tensors and then perform TDE for each sub-tensor. Additionally, we present results for the Tensor-based Eigenfilter. Differently from the CPD-based methods, the Tensorbased Eigenfilter does not require an estimate of the model order to perform TDE of the LOS component. Therefore, the Tensor-based Eigenfilter might be a suitable alternative in a dynamic scenario. Finally, the CPD-GEVD and the HOSVD SECSI methods are not suitable to be combined with the approach in which we create sub-tensors to perform TDE. Since both methods use the dimension of epochs, e.g., dimension K , of the tensors to perform factor matrix estimation in some scenarios withL d > K the tensor factorization becomes impossible. Since the HOSVD SECSI method is not suitable for the sub-tensor approach, we exploit the tensor dimensions that could be used to perform tensor factorization. Thus, we, alternatively, use the proposed Mode 1 HOSVD SECSI with the left-hand slices of the core tensor. In contrast to the HOSVD SECSI, the proposed Mode 1 HOSVD SECSI with left-hand slices of the core tensor uses the antenna dimension, i.e., the third dimension of the tensor Y (d) , to perform tensor factorization while the HOSVD SECSI uses the epochs dimension, i.e., the first dimension of the tensor Y (d) . In Figure 6 we show simulation results for the proposed Mode 1 HOSVD SECSI with left-hand matrix slices with K = 30. As previously described, we use the matrix-based MOS methods to estimate the model order for each epoch, then we group the epochs with the same model order and create sub-tensors. Therefore, since the sub-tensor approach is a solution that aims to create pseudo-static scenarios, we show that the HOSVD SECSI variant is suitable to perform tensor factorization and TDE when applying the sub-tensor approach.

E. TDE CONSIDERING A DYNAMIC SCENARIO WITH MuDe
In this section, we discuss TDE considering a dynamic scenario, a perfectly aligned array of antennas (no errors in the array response model), and the MuDe technique. When applying matrix-based MOS methods, we again grouped the epochs with the same estimated model order, as described previously. The MuDe technique cannot be applied to the Tensor-based Eigenfilter and a simpler matrix-based eigenfilter. For simulations with φ = 10 • most MOS techniques failed to correctly estimate L In Figure 7 we show simulation results for the proposed Mode 1 HOSVD SECSI with left-hand matrix slices applying the MuDe method with K = 30. Note, that the ideal case with known A and is a reference to the smallest error in noisy scenarios. However, in practice, A and must be estimated. Observe that, by inspecting the ideal cases with and without MuDe, if the tensor had higher dimensions, then the gain of denoising would be even higher. We can observe a performance gain when applying MuDe to the sub-tensors, since the MOS methods EFT+FBA+SPS, M-EFT+FBA+SPS, and MDL+FBA+SPS curves show a lower error when combined with MuDe. Besides, since the sub-tensor approach attempts to generate pseudo-static scenarios from dynamic ones, the proposed Mode 1 HOSVD SECSI with the left-hand matrix method achieves an accurate matrix separation, thus achieving improved estimates of the factor matrices.

VI. CONCLUSION
State-of-the-art tensor-based TDE methods are not suitable for scenarios with a time-varying number of multipath components. To overcome this limitation, we proposed a tensor-based framework capable of performing model order estimation and factor decomposition for time-delay estimation in dynamic multipath scenarios. The proposed approach is applicable in time-varying multipath environments and selects the most suitable MOS scheme. To perform high accuracy factor decomposition, we exploit the model order estimates for each slice to group the slices into sub-tensors.
Since, by creating sub-tensors, we achieve pseudo-static sub-scenarios, we obtain a model order larger than the number of epochs. Therefore, the previous state-of-the-art CPD-GEVD and HOSVD SECSI are no longer suitable for dynamic multipath scenarios. We have shown that the proposed Mode 1 HOSVD SECSI provides reliable and accurate estimates when combined with the created sub-tensors. Furthermore, due to the low power of the NLOS components, the NLOS components may be identified as noise when performing MOS. Consequently, the PoD is smaller in case the signals are strongly correlated. Therefore, we have shown that by combining accurate and robust MOS methods with tensor factorization techniques, we obtain more accurate TDE of the LOS signal.
Combining the proposed Mode 1 HOSVD SECSI with MuDe, an additional performance gain can be achieved due to the successful denoising of the generated sub-tensors. Finally, if we increased the number of dimensions and the size of those dimensions, we would obtain an even higher performance gain.