Reconstructing Spatial Localization Error Maps via Physics-Informed Tensor Completion for Passive Sensor Systems

TL;DR

This paper introduces a new method to map localization errors in sensor systems using physics-informed tensor completion, improving accuracy over existing techniques.

Contribution

A novel physics-informed tensor completion framework for reconstructing high-fidelity localization error maps from sparse data.

Findings

The proposed framework improves localization error map accuracy by at least 27.96% over existing methods.

Physics-informed regularization enhances robustness in recovering error maps from incomplete data.

Abstract

Accurate mapping of localization error distribution is essential for assessing passive sensor systems and guiding sensor placement. However, conventional analytical methods like the Geometrical Dilution of Precision (GDOP) rely on idealized error models, failing to capture the complex, heterogeneous error distributions typical of real-world environments. To overcome this challenge, we propose a novel data-driven framework that reconstructs high-fidelity localization error maps from sparse observations in TDOA-based systems. Specifically, we model the error distribution as a tensor and formulate the reconstruction as a tensor completion problem. A key innovation is our physics-informed regularization strategy, which incorporates prior knowledge from the analytical error covariance matrix into the tensor factorization process. This allows for robust recovery of the complete error map even from highly incomplete data. Experiments on a real-world dataset validate the superiority of our approach, showing an accuracy improvement of at least 27.96% over state-of-the-art methods.

Figures

Full text

1. Introduction

Passive localization is ubiquitous in fields ranging from wireless sensor networks (WSNs) to autonomous driving [1,2,3]. The reliability of these systems hinges on a thorough understanding of their performance across the operational domain. Consequently, a high-fidelity spatial localization error map is an indispensable tool, providing crucial insights for system evaluation, operational planning, and strategic sensor deployment. In practice, however, acquiring a complete and dense error map is often infeasible due to cost and logistical constraints, typically yielding only sparse measurements at discrete locations. Therefore, robustly reconstructing the full error distribution from limited data remains a significant challenge in the field.

Currently, the Geometric Dilution of Precision (GDOP) is the most widely adopted metric for performance evaluation [4,5]. The underlying principle is that measurement errors (e.g., in Time-Difference-of-Arrival, TDOA) are amplified by the geometric arrangement of the sensors, where GDOP provides a scalar metric to quantify this amplification factor [6,7]. Typically, unknown model parameters, such as measurement error variances, are estimated statistically from limited emitters and then extrapolated to unobserved regions. To enhance model fidelity in complex environments, numerous extensions to the basic GDOP model have been proposed. These include Weighted GDOP (WGDOP) variants accounting for heterogeneous noise in satellite navigation [8,9,10], as well as adaptations for indoor spaces [11,12], urban canyons [13], UAV-assisted networks [1,14,15], and WSNs [5,16].

Despite these refinements, model-driven methods share a fundamental limitation: they are predicated on idealized assumptions and oversimplified error models. They often fail to capture the complex coupling of error sources in real-world environment ssuch as non-line-of-sight (NLOS) propagation, multipath fading, and clock drift, leading to an unavoidable mismatch between predicted and actual system performance.

To mitigate this model mismatch, data-driven paradigms have emerged as a promising alternative. By learning directly from measurement data, these methods capture fine-grained environmental effects and complex error structures without relying on rigid physical models. Existing data-driven strategies for spatial map reconstruction generally fall into three categories. The first includes classical spatial interpolation techniques, such as Kriging [17,18] and kernel-based methods [19,20,21,22], which estimate values at unobserved locations based on spatial correlation. The second category comprises low-rank matrix and tensor completion methods, including compressed sensing [23], singular value thresholding [24,25,26], and tensor decomposition algorithms [27,28,29]. For instance, Zhang et al. [28] utilized block term decomposition (BTD) to reconstruct electromagnetic maps, while Sun et al. [30] improved accuracy by combining local interpolation with nuclear norm minimization (NNM-T). The third category encompasses machine learning approaches, such as RadioUNet [31], autoencoders [32], and Vision Transformer (ViT)-based methods [33]. While powerful, the efficacy of these purely data-driven methods is critically contingent on the availability of dense and uniformly distributed data [34,35]. In practical scenarios with sparse measurements, their performance degrades sharply, often failing to converge to physically meaningful solutions.

To overcome the limitations of inaccurate physical models and insufficient observational data, this paper proposes a framework that integrates model-based insights with the flexibility of data-driven approaches. First, we model the three-dimensional (3-D) spatial distribution of localization error as a third-order tensor, termed Tensorized GDOP (TGDOP). This representation offers a powerful mathematical tool to describe complex spatial error distributions, capturing intrinsic multi-dimensional structures and inherent anisotropy [36] that scalar metrics like GDOP fail to express. Subsequently, we formulate the reconstruction of the complete error map as a tensor completion problem. To solve this ill-posed problem under sparse measurement conditions, we develop a physics-informed regularization strategy. Specifically, we incorporate prior knowledge from the analytical error covariance matrix directly into the tensor factorization process by imposing polynomial constraints on the factor matrices. This physics-based constraint guides the reconstruction, ensuring the solution adheres to the underlying geometric principles of localization. It effectively compensates for missing observations, dramatically enhancing reconstruction accuracy and robustness without relying on simplified environmental models. Furthermore, this approach is not scenario-dependent and holds potential for extension to complex multi-system fusion scenarios.

Notably, while the concept of physics-informed learning is prominent in Physics-Informed Neural Networks (PINNs) [37], our usage differs fundamentally. PINNs typically integrate partial differential equations (PDEs) into a loss function and require training on large datasets to learn latent physical laws. In contrast, the proposed approach embeds these laws directly into the tensor structure as constraints for decomposition. Consequently, our method operates in a single-shot manner, which is training-free and requires no historical data. Experimental results on a real-world dataset demonstrate that, even with only 1% of observation data available, the proposed framework improves reconstruction accuracy by at least 27.96% compared to state-of-the-art baselines.

The main contributions of this paper are summarized as follows:

• Tensor-Based Spatial Error Modeling (TGDOP): We propose a novel framework that models the spatial distribution of positioning errors as a third-order tensor. Unlike conventional scalar GDOP metrics, this tensor representation explicitly captures the anisotropic characteristics and complex coupling of error sources in real-world 3-D environments.
• Physics-Informed Sparse Reconstruction Algorithm: We develop a robust tensor completion algorithm tailored for extremely sparse observational data. By deriving spatial properties from the theoretical error covariance matrix, we introduce polynomial constraints to the factor matrices during tensor decomposition.
• Training-Free and Model-Robust Performance: The proposed method operates as a single-shot, data-driven approach that does not require historical training data or idealized channel assumptions, validated on both simulated and real-world datasets.

The remainder of this paper is organized as follows: Section 2 describes the problem formulation, and Section 3 provides the preliminaries. The definition and properties of TGDOP are detailed in Section 4, followed by the sparse reconstruction algorithm in Section 5. Section 6 presents simulation and real-data experiments, and conclusions are drawn in Section 7.

2. Problem Statement

Consider a 3-D space discretized into a grid of $[eqn]$ cells. The coordinate corresponding to a grid cell with index $[eqn]$ is denoted by $[eqn]$ , where $[eqn]$ represents the set of all such indices, with  cardinality $[eqn]$ . Let $[eqn]$ be the set of indices for grid cells containing emitters.

The characterization of positioning error begins with understanding its statistical nature at a specific location $[eqn]$ . Assume there are N independent positioning results $[eqn]$ (for $[eqn]$ ) for an emitter at a true location $[eqn]$ , then the positioning error for each result is given by $[eqn]$ . The covariance matrix of the positioning error at $[eqn]$ is then calculated by

$$[eqn]$$

where $[eqn]$ . Suppose the positioning error measurements are available for any grid cell, then a heatmap of the spatial distribution of localization error can be generated by calculating the Root Mean Square Error (RMSE) value at each grid point, $[eqn]$ , as  depicted in Figure 1.

In real-world applications, however, the acquisition of valid measurements is restricted to a very limited number of locations equipped with emitters. This inherent limitation means that statistical approaches can only characterize the positioning error locally and are incapable of perceiving the error map across the entire space. To overcome this, conventional methods use a parameterized empirical model for $[eqn]$ instead of a statistical one in (1), allowing extrapolation to emitter-free locations. These models typically define a few unknown parameters based on key error sources, which are then solved for using the limited available measurements. A typical example is the error map derived by the GDOP model for TDOA systems, which primarily accounts for errors in TDOA measurement noise $[eqn]$ and sensor position uncertainties $[eqn]$ :

$$[eqn]$$

Here, the specific forms of the geometry matrix $[eqn]$ and the intermediate parameter error covariance $[eqn]$ are detailed in Section 4.2. The variance of the l-th TDOA measurement, $[eqn]$ , can be further expressed based on signal parameters [38] as

$$[eqn]$$

where $[eqn]$ is the root mean square (RMS) signal bandwidth, B is the signal bandwidth, T is the RMS integration time, and  $[eqn]$ is the effective signal-to-noise ratio for the l-th TDOA calculation.

The parametric models can align with statistical estimates calculated by (1) if the underlying model assumptions are accurate. However, real-world scenarios involve complex propagation channels and diverse, often unmodeled, error sources. For instance, in typical NLOS scenarios, while the factors detailed (such as signal bandwidth and SNR) in (3) influence the standard deviation of TDOA measurements, the error component arising from NLOS propagation often plays a more dominant and decisive role. However, the GDOP model would be highly inaccurate because it does not account for NLOS errors. Moreover, the  need to design distinct models for $[eqn]$ tailored to different localization systems and specific application scenarios significantly curtails the practical utility and generalizability of conventional models.

To address the limitations of conventional models, the tensor-based framework surpasses conventional models by learning complex spatial distributions directly from data, independent of empirical models. We define the tensor representation of the 3-D positioning error map as $[eqn]$ . Its components along the Cartesian axes (x, y, and z) are denoted as $[eqn]$ , $[eqn]$ , and $[eqn]$ , respectively. Then, the relationship between $[eqn]$ and its components is given by

$$[eqn]$$

In practical scenarios with spatially sparse emitters, the error statistics can only be measured at indices $[eqn]$ . The measurement model for an element at index $[eqn]$ is

$$[eqn]$$

where $[eqn]$ and $[eqn]$ represent the observation errors associated with $[eqn]$ and $[eqn]$ , respectively.

The learnable parameters in traditional models are limited to a predefined set of error sources from a prior model. In contrast, the tensor model’s parameters consist of the factor matrices and the core tensor (detailed in Section 4.1), where the factor matrices characterize the axial trends of the error distribution, and the core tensor represents the coupling relationships among these axial components. Thus, the tensor model can better adapt to complex real-world environments for inference in regions with no observations.

In this context, our primary problem is to reconstruct the complete underlying error component tensors $[eqn]$ (and, consequently, the total error tensor $[eqn]$ ) across the entire spatial domain $[eqn]$ from sparse and noisy measurements $[eqn]$ (and $[eqn]$ ) available only at indices $[eqn]$ . It is typically the case that the number of observed points is much smaller than the total number of grid points, i.e.,  $[eqn]$ , and the observed locations $[eqn]$ are often randomly distributed.

Notation: In this paper, we follow the established convention in signal processing. We denote scalars, vectors, matrices, and tensors with lowercase letters $[eqn]$ , boldface lowercase $[eqn]$ , boldface capitals $[eqn]$ , and underlined boldface capitals $[eqn]$ , respectively. $[eqn]$ denotes the $[eqn]$ -th element of a third-order tensor. The colon notation indicates the sub-tensors of a given tensor. Italic capitals are also used to denote the index upper bounds. We use the superscripts $[eqn]$ and $[eqn]$ to represent the inverse and the pseudo-inverse of a matrix, respectively. In addition, we use $[eqn]$ to denote the Frobenius norms of a tensor. The operator $[eqn]$ stacks its scalar arguments in a square diagonal matrix, $[eqn]$ is the covariance of $[eqn]$ and $[eqn]$ , and  $[eqn]$ computes the trace of a matrix. The  $[eqn]$ operator is defined as $[eqn]$ . The outer product follows that $[eqn]$ . The Hadamard product follows that $[eqn]$ . The Kronecker product is derived by $[eqn]$ , and its capital symbol $[eqn]$ is used to denote the multiplicative form. Then, the Khatri–Rao product is defined as the block-wise Kronecker product: $[eqn]$ . The $[eqn]$ denotes the Kronecker delta. The  $[eqn]$ identity matrix is represented by $[eqn]$ , $[eqn]$ zero matrix is denoted by $[eqn]$ , and  $[eqn]$ is a column vector of all ones of length N.

3. Preliminaries

We begin with fundamental tensor analysis concepts to grasp the approach outlined in this paper. The issue discussed in this paper aligns with BTD, so the emphasis is on introducing the key concepts of BTD.

3.1. Mode-N Unfolding and Mode Product

For an N-th order tensor $[eqn]$ , the mode-n unfolding of $[eqn]$ , denoted as $[eqn]$ , unfolds the n-th fiber as the columns of the resulting matrix. It maps the tensor element $[eqn]$ to the matrix element $[eqn]$ , where

$$[eqn]$$

The mode-n product of the $[eqn]$ with the $[eqn]$ is denoted as $[eqn]$ , and it follows that

$$[eqn]$$

Its matrix representation is expressed as

$$[eqn]$$

3.2. Block Tensor Decomposition in Multilinear Rank-(L,M,N) Terms

Unlike matrices, tensor rank determination is an NP-hard problem, leading to varied rank definitions among algorithms. The concept of multilinear rank is proposed in BTD. The tensor’s mode-n rank is defined as

Definition 1(Mode-n rank). The mode-n rank of a tensor $[eqn]$ is the dimension of the subspace spanned by its mode-n vectors.

Then, a third-order tensor’s multilinear rank is $[eqn]$ if its mode-1 rank, mode-2 rank, and mode-3 rank are equal to L, M, and N, respectively.

A decomposition of a tensor $[eqn]$ in a sum of rank-(L,M,N) terms can be written as

$$[eqn]$$

where $[eqn]$ is rank-(L,M,N), and $[eqn]$ , $[eqn]$ , and  $[eqn]$ are matrices with full column rank. Typically, the symbol $[eqn]$ is referred to as the core tensor, while the matrices $[eqn]$ , $[eqn]$ , and  $[eqn]$ are denoted as the factor matrices. Calculating the mode-n unfolding of (9), we have

$$[eqn]$$

where $[eqn]$ , $[eqn]$ , and  $[eqn]$ .

4. TGDOP and Its Properties

In this section, we introduce the TGDOP to characterize the distribution of localization errors and investigate the potential properties of TGDOP by analyzing the covariance matrix of positioning errors. Initially, we present a decomposed formulation of TGDOP, which explicitly accounts for the anisotropic characteristics inherent in localization errors. Subsequently, we analyze the general expressions for the error covariance matrix under various localization schemes. This analysis aims to establish the universality of the TGDOP model across different localization methodologies. Finally, by mapping the properties of the covariance matrix to the fundamental characteristics of the factor matrices within the TGDOP tensor space, we discuss the theoretical underpinnings that enable data-driven reconstruction using this model.

4.1. Tensor Model of Positioning Error Distribution

In real-world applications, conventional positioning error models are disturbed by model mismatch because their constrained parametric form limits their ability to capture the multiplicity of error sources found in complex environments. In contrast, tensor models are structured with far greater expressive power, allowing them to be effectively data-driven in constructing models that are sufficiently accurate and adaptable for practical applications.

The expressive power of tensor models stems from their core tensor and factor matrices. Analogous to the singular value decomposition (SVD) in 2-D cases, the 3-D tensor $[eqn]$ $[eqn]$ could be characterized by eigenspaces defined by three sets of eigenvectors. We denote these eigenspaces as $[eqn]$ , $[eqn]$ , and $[eqn]$ , corresponding to the factor matrices on the x, y, and z dimensions, respectively. Consequently, we have

$$[eqn]$$

where $[eqn]$ is the core tensor of $[eqn]$ . Substituting (11) into (4), $[eqn]$ can be expressed as

$$[eqn]$$

The theoretical TGDOP model is presented in (12). Employing a BTD framework, this model effectively maps the spatial distribution of the diagonal components of the covariance matrix $[eqn]$ onto a set of latent core tensors, $[eqn]$ . When this formulation is integrated with the expression for TGDOP measurements given in (5), it yields $[eqn]$ . The probability distribution of the measurement noise $[eqn]$ here satisfies Theorem 1.

Theorem 1(Measurement Noise Distribution). Consider a 3-D scenario where an emitter is located at $[eqn]$ , and there are N positioning results for this emitter. Let $[eqn]$ represent the covariance matrix of the positioning error. When N is relatively larger than the diagonal elements in $[eqn]$ , the observation noise of $[eqn]$ and $[eqn]$ at $[eqn]$ follows that

$$[eqn]$$

The proof is relegated to Appendix A. Theorem 1 establishes the condition under which measurements conform to a Gaussian distribution.

Notably, the TGDOP model presented in (12) involves not only the tensor $[eqn]$ but also all its constituent core tensors and factor matrices. These elements enable the reconstruction of both $[eqn]$ and its components, $[eqn]$ . In essence, TGDOP can be represented by a higher-order tensor $[eqn]$ satisfying $[eqn]$ . We demonstrate that $[eqn]$ , which lacks directional information, is a collapse of the directionally informed $[eqn]$ along different dimensions. The relationship between $[eqn]$ and $[eqn]$ is further elucidated in Appendix C through the definition of a tensor operation.

4.2. The Covariance Matrix of Positioning Error

The covariance matrix of the positioning error is directly related to the TGDOP measurements, so the properties of the covariance matrix are the theoretical basis for designing the TGDOP factor matrix. Many documents have derived the covariance matrices of positioning errors under different positioning systems [36,39,40,41,42] in 2-D cases. This subsection proposes the covariance matrices of positioning errors for the TDOA positioning system in a 3-D Cartesian coordinate system.

Now, by defining the covariance matrix of the positioning error as $[eqn]$ , we have $[eqn]$ for a 3-D scenario. Assume there are N sensors for TDOA positioning, where the main sensor and the auxiliary sensors are located at $[eqn]$ and $[eqn]$ , respectively, and the radiation source is located at $[eqn]$ . Let $[eqn]$ denote the TDOA measurements for the i-th group of sensors, $[eqn]$ denote the variance of $[eqn]$ , and  $[eqn]$ denote the variance of the sensor position error. The i-th row of the direction cosine matrix is given by

$$[eqn]$$

Let c denote the light speed, then the covariance matrix of the positioning error follows that $[eqn]$ , where

$$[eqn]$$

Next, by defining $[eqn]$ to represent the sum of the covariance matrices of all possible estimation errors, we have

$$[eqn]$$

It is found that under different localization regimes, the error covariance matrix $[eqn]$ can consistently be derived from $[eqn]$ and $[eqn]$ in a form analogous to (15). This is because (15) elucidates the principal sources of localization error. Specifically, $[eqn]$ represents the inherent estimation errors within the localization system, while $[eqn]$ quantifies the anisotropic amplification effect on these errors attributable to the station geometry. Therefore, (15) can be considered a general expression for the localization error covariance matrix. Moreover, according to (5), the TGDOP measurements also incorporate this universal expressive capability, which explains its applicability across various localization systems.

4.3. Properties of TGDOP Derived from Error Covariance Matrix

A thorough understanding of its factor matrix properties is paramount to enable the effective reconstruction of the TGDOP model introduced in (12). This subsection investigates these properties by first examining the positioning error covariance matrix. We then translate the characteristics of this covariance matrix to the factor matrices in the tensor space, thereby establishing the smoothness, non-negativity, and low-rank properties of the proposed TGDOP model.

4.3.1. Spatial Smoothness

The mathematical notation $[eqn]$ is commonly used to describe the smoothness of a function. Specifically, $[eqn]$ indicates that a function is continuous over its domain, while $[eqn]$ signifies that its first derivative exists and is also continuous over its domain (i.e., the function is continuously differentiable). Theorem 2 outlines the conditions for spatial continuity when the covariance matrix is treated as a matrix function of spatial position. A detailed derivation is available in Appendix B.

Theorem 2 (Continuity of Derivatives).

• $[eqn]$ when $[eqn]$ satisfies the following conditions:*

• • $[eqn]$ ;*
• • $[eqn]$ .*

(x, y, z) represent the position variables of the radiation source in the Cartesian coordinate system.

Theorem 2 illustrates that the spatial continuity of $[eqn]$ plays a crucial role in the TGDOP model. Subsequently, we will analyze the continuity conditions for $[eqn]$ in the TDOA localization system. Calculate the $[eqn]$ :

$$[eqn]$$

where $[eqn]$ is defined in (13). It can be readily observed that for an interval I where the condition $[eqn]$ holds, the partial derivative $[eqn]$ is continuous on I (i.e., $[eqn]$ ). The forms of $[eqn]$ and $[eqn]$ are analogous to that of $[eqn]$ , and thus, they adhere to the same continuity properties. Furthermore, $[eqn]$ itself is continuous on I (i.e., $[eqn]$ ). Additionally, under a typical/standard sensor configuration, $[eqn]$ possesses full column rank. Consequently, it can be readily concluded that $[eqn]$ is based on Theorem 2.

According to (5) and (12), the property that $[eqn]$ can be mapped to the tensor space. This implies that the dimension of the factor matrices’ solution space in (12) can be reduced through structured constraints, which is important for subsequent sparse reconstruction algorithms.

4.3.2. Non-Negative

Theorem 3 states that the covariance matrix $[eqn]$ is symmetric and positive semi-definite. Let $[eqn]$ . Then, for all $[eqn]$ , we have

$$[eqn]$$

Considering the definition of $[eqn]$ for $[eqn]$ in (5), it follows that $[eqn]$ . By defining $[eqn]$ and $[eqn]$ , respectively, similar conclusions can be drawn for $[eqn]$ and $[eqn]$ . Leveraging this property, non-negative constraints are applied to the factor matrices of (12) in the subsequent algorithm. This confines the solution space to the domain of non-negative real numbers, reducing the variable dimensions in parameter optimization.

Theorem 3 (Symmetric Positive Semi-definite). Let $[eqn]$ , where $[eqn]$ denotes the positioning error vector. Then, $[eqn]$ satisfies the following:

• • $[eqn]$ ;*
• • $[eqn]$ .*

Proof. 

$$[eqn]$$

   □

4.3.3. Low Rank

To investigate the hypothesized low-rank property of the proposed tensor model $[eqn]$ , which is constructed from individual covariance matrices $[eqn]$ as per (5), our approach is to analyze the rank characteristics of its constituent slices using SVD. The methodology involves simulating a scenario with high variability in the underlying error parameters to assess if low-rank structures persist even under such complex conditions.

The positioning error covariance matrix $[eqn]$ at a location $[eqn]$ is generally expressed from (15) as $[eqn]$ . Here, $[eqn]$ is a geometry-dependent transformation matrix, and  $[eqn]$ is a symmetric, positive, semi-definite matrix encapsulating the statistics of underlying sensor measurements or intermediate error parameters, influenced by environmental factors $[eqn]$ . Performing a Cholesky decomposition on $[eqn]$ yields $[eqn]$ , where $[eqn]$ is a lower triangular matrix. Consequently, the covariance matrix can be written as

$$[eqn]$$

The matrix $[eqn]$ is determined by the fixed station configuration. Thus, the variability and specific structure of the positioning error map across different environments are primarily determined by $[eqn]$ .

Assuming there are $[eqn]$ independent intermediate variables contributing to $[eqn]$ (e.g., if N is the number of variables per sensor and M is the number of sensors or total measurements, then $[eqn]$ could be $[eqn]$ ), $[eqn]$ is an $[eqn]$ matrix, and thus $[eqn]$ . We populate the non-zero elements of $[eqn]$ with random values to simulate a complex error environment lacking strong prior structural information. The tensor $[eqn]$ is then constructed using these generated $[eqn]$ matrices. We then perform SVD on matrix slices of $[eqn]$ . For instance, if a $[eqn]$ results in a matrix $[eqn]$ (for $[eqn]$ ), the normalized cumulative energy is defined as $[eqn]$ , where $[eqn]$ is the j-th singular value of $[eqn]$ in descending order, and  $[eqn]$ . Then, the average of cumulative energy for K slices is defined as $[eqn]$ .

Figure 2 plots $[eqn]$ versus i for such an analysis, where the considered tensor slice dimensions result in $[eqn]$ singular values. It is found from Figure 2 that the first five singular values account for the majority of the energy in the analyzed tensor. This rapid decay of singular values suggests that the matrix representation of the positioning error distribution, even under these simulated complex conditions, exhibits a strong low-rank characteristic. This finding provides empirical support for the premise that the spatial distribution of positioning error can be effectively modeled by low-rank structures, offering a foundation for the low-rank reconstruction algorithm developed in this paper. While this analysis is based on a simulated scenario, the results suggest a general tendency that is further exploited by our proposed method.

5. Reconstruction Algorithm for TGDOP

In this section, we propose the sparse reconstruction algorithm for TGDOP. In cases where the positioning error sources are unidentified and limited emitters are available, we propose a BTD-based algorithm to enhance reconstruction performance under sparse measurements by constraining the factor matrix. In addition, we extend the algorithm when the directional measurements of positioning error are available. Specifically, we constrain the corresponding factor matrices with the directional measurements of positioning error in the objective function and reconstruct the optimized factor matrices into TGDOP.

5.1. BTD-Based Sparse Reconstruction Algorithm for TGDOP

We begin by outlining the context for the proposed algorithm. In this paper, a “sparse observational scenario” signifies a situation where calibration emitters are present in only a limited number of grids within the target area. Sensors receive electromagnetic signals from these emitters. Subsequent statistical analysis of multiple positioning results at each grid yields positioning errors; the squared magnitudes of these errors constitute the observational values, as detailed in (5).

As outlined in (12), the primary objective is the reconstruction of the complete positioning error tensor. This goal is achieved by solving the following optimization problem, which models the error tensor as a sum of R BTD components:

$$[eqn]$$

In this formulation, $[eqn]$ denotes the measurement tensor. Each term in the summation, $[eqn]$ , represents a BTD component comprising a core tensor $[eqn]$ and associated factor matrices $[eqn]$ , $[eqn]$ , and $[eqn]$ . In accordance with (12), the number of BTD blocks is set to 3, representing the three directional components: x, y, and z. The dimensions $[eqn]$ satisfy $[eqn]$ (for $[eqn]$ ), and these factor matrices are constrained to possess full column rank, a common requirement for model identifiability.

5.1.1. Designing the Structure of Factor Matrices

Section 4.3 analyzes the intrinsic properties of the positioning error tensor $[eqn]$ . The factor matrices in the tensor decomposition represent projections of $[eqn]$ ’s features along its different modes. Consequently, the  characteristics of $[eqn]$ should inform the structure of these factor matrices. Theorem 2 reveals the spatial smoothness inherent in positioning errors. This implies that the column vectors of the factor matrices (which capture variations along spatial dimensions) can effectively represent most of the information using low-order polynomials. This prior knowledge is particularly valuable in scenarios with sparse observational data. Standard random initialization of factor matrices fails to incorporate this fundamental prior, potentially leading to models that are poorly constrained by sparse measurements and thus increasing the risk of convergence to suboptimal local minima.

Theorem 2 establishes the continuity of elements within the tensor of positioning error distribution. Here, we relate this property to the factor matrices of the tensor via Theorem 4, the detailed proof of which is provided in Appendix E.

Theorem 4(Vandermonde-like structure). Given that the error distribution tensor $[eqn]$ exhibits $[eqn]$ continuity with respect to spatial position and possesses a low-rank structure, it follows that the column vectors of its factor matrices reside in a subspace generated by low-order polynomials.

Based on Theorem 4, we propose structuring the column vectors of each factor matrix using polynomial generators. The optimization then focuses on two main steps. First, a  polynomial coefficient matrix is randomly initialized. Second, these coefficients are used to construct the columns of the actual factor matrix. Subsequent optimization iterations refine these polynomial coefficients rather than directly adjusting the elements of the factor matrices. This methodology explicitly incorporates the prior knowledge of spatial smoothness, promoting spatial continuity in the reconstructed tensor. Furthermore, by embedding this structural prior, the approach reduces the search space and guides the optimization, thereby decreasing the likelihood of converging to undesirable local optima.

Consider optimizing $[eqn]$ in (19) as an illustrative example. We first randomly initialize a polynomial coefficient matrix $[eqn]$ , where m is the number of coefficients used to define a polynomial of degree $[eqn]$ . Each of the $[eqn]$ columns of $[eqn]$ is generated using a distinct set of m polynomial coefficients. The initial factor matrix $[eqn]$ is then constructed as

$$[eqn]$$

where $[eqn]$ is the Vandermonde matrix whose columns form a basis for polynomials up to degree $[eqn]$ . Since $[eqn]$ is full column rank (given $[eqn]$ ), optimizing the objective function with respect to $[eqn]$ is equivalent to optimizing with respect to $[eqn]$ . The optimization task thus shifts from determining $[eqn]$ (with $[eqn]$ parameters) to determining $[eqn]$ (with $[eqn]$ parameters), significantly reducing the dimension of parameters to be optimized in the solution space when $[eqn]$ . This reduction effectively smooths the optimization landscape and prevents the algorithm from getting trapped in bad local minima associated with high-frequency noise, thereby improving the robustness to random initialization.

To ensure the non-negativity of the factor matrices indicated in Theorem 3, we define the final factor matrix $[eqn]$ through an element-wise squaring operation:

$$[eqn]$$

Upon substituting Equations (20) and (21) into the objective function (denoted as $[eqn]$ ) defined in (19), the optimization problem concerning the factor matrix $[eqn]$ can be re-expressed as

$$[eqn]$$

Consequently, optimizing the objective function $[eqn]$ from (19) concerning the factor matrix $[eqn]$ (now denoted $[eqn]$ to reflect its construction) becomes equivalent to optimizing for the polynomial coefficient matrix $[eqn]$ . For brevity, any subsequent reference to optimizing a factor matrix implies the optimization of its underlying polynomial parameter matrix $[eqn]$ .

5.1.2. Sparse Formulation of the Objective Function

In practical applications, tensor reconstruction often addresses scenarios with sparse measurements, where the measured data $[eqn]$ covers only a minor portion of the entire target area. To model this, we introduce a binary sampling tensor $[eqn]$ , where $[eqn]$ if the element $[eqn]$ is observed (i.e., $[eqn]$ ), and 0 otherwise. Assuming the unobserved entries in $[eqn]$ are represented as zeros, the optimization problem for sparse data is formulated as

$$[eqn]$$

5.1.3. Solving Equation (23) Using Block Coordinate Descent

Now, we employ a Block Coordinate Descent (BCD) strategy to solve (23), which involves iteratively minimizing the objective function with respect to one block of variables while keeping the others fixed. This leads to the following four sub-problems, derived from the mode-n matricized forms of the objective function:

$$[eqn]$$

where

$$[eqn]$$

The variables to be optimized in (23) are the sets of factor matrices $[eqn]$ , $[eqn]$ , $[eqn]$ (collectively denoted $[eqn]$ ), and the set of core tensors $[eqn]$ . Denoting the objective function in (23) by $[eqn]$ , the  BCD method iteratively solves the following:

$$[eqn]$$

Each sub-problem in (26) is a linear least squares problem concerning the optimized variable block and is therefore convex. Consequently, an Alternating Least Squares (ALS) approach, a specific type of BCD, can be employed, as detailed in Algorithm 1. Algorithm 1 ALS algorithm for solving (23).

• Initialize $[eqn]$ .
• while not converged do
• for  $[eqn]$  do
• Update  $[eqn]$ :
• $[eqn]$ .
• for  $[eqn]$  do
• $[eqn]$ .
• $[eqn]$ .
• end for
• for  $[eqn]$  do
• Perform QR-factorization: $[eqn]$ .
• $[eqn]$ , $[eqn]$ .
• end for
• end for
• Update  $[eqn]$ :
• $[eqn]$ , $[eqn]$ .
• $[eqn]$ .
• end while
• Reconstruct $[eqn]$ and $[eqn]$ using converged factors via (11) and (12), respectively.
• Return $[eqn]$ and $[eqn]$ .

The focus of this study lies in tensor modeling methodology, so we employed a straightforward ALS-BCD optimization scheme without an exhaustive comparison of alternative reconstruction strategies. Figure 3 illustrates the average convergence curves of the proposed algorithm and the classical BTD algorithm in LOS and NLOS scenarios, based on 100 random initializations. In the figure, the x-axis represents the number of optimization steps performed by the ALS algorithm, and the y-axis represents the total squared Frobenius norm error defined in (23). This value is calculated based on the normalized measurement tensor and is used to illustrate the relative magnitude of the convergence trend. It can be observed that the proposed algorithm outperforms the classical BTD method in terms of both convergence speed and stability, and exhibits sufficient robustness to initial value selection. Nevertheless, to extend the algorithm to larger tensor spaces, advanced techniques such as momentum-based updates, stochastic optimization, or distributed computation can be incorporated to further boost convergence speed and robustness.

For the aforementioned solution process, the derivation of its computational complexity is provided in Appendix D:

$$[eqn]$$

where $[eqn]$ denotes the number of observed samples. It is evident that the complexity of the proposed algorithm is acceptable under the low-rank assumption. However, when the multilinear rank is large, reconstruction accuracy must be sacrificed to reduce time costs. Future research may overcome this technical bottleneck by modifying the optimization strategy.

5.2. Solution Approach with Available Directional Measurements

The rotational ambiguity of sub-tensors is a known challenge in tensor decomposition and is often handled with prior constraints, initialization, or post-processing [28]. This subsection explores an alternative solution strategy when directional components of the positioning error are available as measurements. It is emphasized that the proposed method does not decompose the tensor into three arbitrary blocks. Instead, it reconstructs three physically meaningful components grounded in the definition of the error covariance matrix. This ensures that the physical correspondence is guaranteed by the model’s design from the outset, rather than being an ambiguous interpretation of the decomposition results.

While existing research on positioning errors has largely concentrated on scalar distributions (e.g., GDOP-based models) that overlook directional characteristics, such information is not irretrievable. When a set of N positioning measurements $[eqn]$ is collected for an emitter at a known location $[eqn]$ , the directional error components $[eqn]$ can be extracted through statistical analysis.

The availability of these directional measurements $[eqn]$ enables a refined modeling approach. We can reformulate the objective function to model each of the R observed directional tensors with a dedicated BTD component, leading to the expression:

$$[eqn]$$

This objective would then be decomposed into sub-problems:

$$[eqn]$$

Utilizing directional measurements aims to provide richer constraints for the factor matrices and core tensors, potentially improving the accuracy of directional positioning error reconstruction. Similar to Algorithm 1, a BCD approach can be applied to solve these convex sub-problems; the detailed algorithmic steps are omitted here for brevity.

5.3. Lower Bound on Sample Complexity

In this section, we discuss the lower bound of the data requirement of the proposed algorithm under the random uniform sampling strategy. We consider the recovery of a K-th order tensor $[eqn]$ , which is represented as a sum of R low-rank Tucker tensors $[eqn]$ : $[eqn]$ . Each component tensor $[eqn]$ has a Tucker decomposition:

$$[eqn]$$

where $[eqn]$ is the core tensor and $[eqn]$ is the factor matrix for component r along mode k, with $[eqn]$ being the mode-k Tucker rank. The analysis results on Tucker decomposition [43,44] indicate that the lower bound for the sample complexity is determined by the most difficult-to-recover matrix unfolding, $[eqn]$ . The recovery difficulty of a matrix is lower-bounded by a function proportional to $[eqn]$ , that is,

$$[eqn]$$

In Section 5.1.1, we assume that each column of every factor matrix $[eqn]$ is a polynomial sampled at $[eqn]$ distinct points. This structural prior can be formulated as a factorization $[eqn]$ , where $[eqn]$ and $[eqn]$ . Here, $[eqn]$ represents the mode-k rank of the r-th component of the tensor, and $[eqn]$ represents the order used for polynomial approximation. Each row of $[eqn]$ contains the coefficients for the corresponding column of $[eqn]$ in the basis defined by $[eqn]$ .

The mode-k unfolding of the r-th component tensor, $[eqn]$ , is given by

$$[eqn]$$

Substituting the polynomial constraint $[eqn]$ for all modes $[eqn]$ , we obtain

$$[eqn]$$

The expression for $[eqn]$ shows that its column space is a subspace of the column space of the known matrix $[eqn]$ :

$$[eqn]$$

The search for the left singular vectors is thus restricted from the ambient $[eqn]$ space to the $[eqn]$ -dimensional subspace spanned by the columns of $[eqn]$ . The effective row dimension is reduced from $[eqn]$ to $[eqn]$ . Using the mixed-product property of the Kronecker product on the right-hand side of the expression:

$$[eqn]$$

This implies that the row space of $[eqn]$ is contained within the column space of the matrix $[eqn]$ , whose dimension is $[eqn]$ . The effective column dimension is therefore reduced from $[eqn]$ to $[eqn]$ .

For the matrix recovery problem of $[eqn]$ , we conclude that it is rank- $[eqn]$ with effective $[eqn]$ row dimension and $[eqn]$ column dimension. Applying the standard lower bound for matrix recovery yields the bound for this specific unfolding:

$$[eqn]$$

Theorem 5 gives the overall sample complexity lower bound for the BTD recovery problem, showing that the sample complexity is determined by the bottleneck among all $[eqn]$ unfoldings.

Theorem 5. For the recovery of a BTD tensor with polynomial factor matrix constraints, the required number of measurements m has the following lower bound:

$$[eqn]$$

where $[eqn]$ is the Tucker rank and $[eqn]$ is the degree of freedom for the factor matrix of component r along mode k.

Figure 4 shows the log-reconstruction error under different observation ratios. The x-axis represents a ratio between 0 and 1, indicating the proportion of observed data points relative to the total number of tensor elements. The y-axis denotes the size N of the $[eqn]$ three-dimensional tensor. The colorbar of the heatmap represents the reconstruction error on a logarithmic scale. The performance of the conventional BTD algorithm and the proposed method is presented in Figure 4a and Figure 4b, respectively. The results indicate that under identical sampling conditions, our method exhibits superior performance by achieving a lower reconstruction error. This consequently reduces the sample lower bound for achieving high-fidelity tensor reconstruction.

6. Results with Measurements and Simulations

To evaluate the performance and robustness of the proposed algorithm under diverse conditions, we conduct simulations on two typical TDOA-based positioning scenarios, denoted as $[eqn]$ and $[eqn]$ , respectively. The experimental findings and conclusions can also be extended to Direction-of-Arrival (DOA) positioning scenarios. Nevertheless, specific results from DOA experiments are not further detailed in this paper for brevity.

The scenario $[eqn]$ utilizes a 3-D TDOA positioning system with four sensors. Three sets of TDOA measurements are the intermediate parameters. The standard deviations for these TDOA measurement errors are set to 18 ns, 20 ns, and 25 ns. The standard deviation of sensor position errors is modeled to be 0.5 m. The sensor coordinates (in meters) are $[eqn]$ , $[eqn]$ , $[eqn]$ , and $[eqn]$ , forming a “Y”-shaped constellation. This configuration is chosen not only because it represents the most common and high-performing geometry in distributed positioning, but also to align with the experimental setup described later. Due to a shared reference sensor in TDOA calculations, correlations are introduced among the TDOA measurement errors. Let the TDOA error vector be $[eqn]$ . The correlation matrix $[eqn]$ for these TDOA errors has $[eqn]$ . The maximum simulated detection range extends to 8 km, with the target area’s elevation spanning from 0.5 km to 5 km. This target volume is discretized into a grid of $[eqn]$ points. Scenario $[eqn]$ builds upon $[eqn]$ by incorporating multipath errors within a quarter of the region spanning elevations from 0.5 km to 1.4 km. The multipath error affecting the TDOA estimates is modeled using an exponential distribution with a mean of 100 ns.

To quantitatively assess the algorithm’s performance, several error metrics were employed. Mean Frobenius Norm Error (MFNE) is used to evaluate the average reconstruction error of a tensor $[eqn]$ over multiple trials: $[eqn]$ , where $[eqn]$ is the true tensor, $[eqn]$ is the tensor reconstructed in the m-th Monte Carlo trial, and M is the total number of trials. Relative Frobenius Norm Error (RFNE) measures the reconstruction error relative to the magnitude of the true tensor: $[eqn]$ . This metric is typically calculated for each trial or averaged over M trials. Signal-to-Noise Ratio (SNR) represents the ratio of the true signal power to the noise power in the measurements, defined in decibels (dB) as $[eqn]$ , where $[eqn]$ is the true underlying (noise-free) tensor and $[eqn]$ is the observed noisy tensor.

6.1. Multilinear Rank Analysis

In this subsection, we employ the L-curve method (Note: The L-curve corner identifies the inflection point where further increases in model parameters yield diminishing returns in error reduction; the rank at this corner is selected as the optimal multilinear rank to balance accuracy and model simplicity.) to determine a suitable multilinear rank for the tensor model optimized by our proposed algorithm. Given that the concept of low rank is relative to a tensor’s dimensions, we estimate the multilinear rank for the target tensor.

To manage computational load during rank analysis, especially for large tensors, we introduce a down-sampling factor $[eqn]$ . The analysis is performed on a down-sampled version of the target data tensor, denoted as $[eqn]$ . The L-curve method typically involves plotting a measure of reconstruction error against model complexity. The corner of this curve often indicates an optimal trade-off. We identify this region by conducting trials with various combinations of multilinear ranks and select the ranks near the L-curve’s inflection point. This selection is further guided by the criterion that the RFNE between the $[eqn]$ and its low-rank approximation $[eqn]$ falls below a predefined threshold $[eqn]$ .

Figure 5 illustrates the estimation of the tensor’s multilinear rank using the L-curve method at $[eqn]$ . The x-axis represents the compression ratio, defined as the percentage of parameters required for reconstruction relative to the total tensor size, while the y-axis denotes the RFNE. The selection of $[eqn]$ provides an optimal trade-off, because it reduces the number of elements by a factor of $[eqn]$ , enabling rapid iterative rank searching while maintaining the macro-scale spatial correlation of the error distribution as guaranteed by the continuity of the TGDOP model. Furthermore, extra experimental results also reveal that across various down-sampling factors $[eqn]$ from 1 to 20, the estimated multilinear rank for each mode consistently remained below 10% of that mode’s original dimension, underscoring the inherent low-rank nature of the data.

6.2. Performance Under Sparse Measurements

This subsection evaluates the algorithm’s reconstruction performance when the emitters are sparsely and randomly distributed within the target area. We define the missing emitter ratio ( $[eqn]$ ) as the proportion of grid points in the target area lacking emitter data [45]. We conducted experiments with $[eqn]$ varying from 80% to 99% (corresponding to an observation ratio of 20% to 1%) to assess performance under increasing sparsity where the observation ratio is defined as 100% MER. For these experiments, the down-sampling factor $[eqn]$ is set to 10, the SNR is maintained at 22 dB, and the number of Monte Carlo trials (M) is 50.

Figure 6 compares the reconstruction performance of the proposed algorithm against several baseline methods across different MER levels. The baseline algorithms included in the comparison encompass both classic reconstruction algorithms and deep learning (DL)-based generative algorithms. It is worth noting that the scenario addressed in this paper is a “single-shot reconstruction” problem, characterized by extremely sparse valid observations and the absence of a historical database. Traditional DL methods typically fail in such zero-shot and extremely sparse conditions due to the inability to acquire sufficient training data. To enable the neural networks to function in this context, we generated 50 sets of historical data to serve as prior information for each training session. Although this constitutes an unfair comparison (as it provides an advantage to the DL baselines), we include these results to verify the superiority of the proposed algorithm from another perspective.

Figure 6a,b demonstrate that the proposed algorithm achieves the lowest $[eqn]$ , indicating superior reconstruction performance. Specifically, at a high sparsity level of $[eqn]$ (only 1% observation ratio), the proposed algorithm exhibits a performance improvement ranging from 43.58% to 97.23% compared to the baseline methods, as detailed in Table 1. Here, $[eqn]$ is defined as $[eqn]$ . Note that the GeneralBTD method represents the standard BTD without the proposed physics-based polynomial constraints. The significant performance gap between the proposed algorithm and GeneralBTD quantifies the contribution of the proposed physical constraints. Furthermore, comparative analysis indicates that the proposed algorithm attains comparable reconstruction accuracy while requiring significantly fewer emitters. This substantial reduction in required emitters can translate to considerable operational cost savings.

Figure 6c,d present boxplots illustrating the distribution of $[eqn]$ at different elevations for $[eqn]$ . As shown in Figure 6d, the introduction of NLOS error at altitudes between 0.5 km and 1.4 km leads to significant reconstruction errors for the WGDOP algorithm in this range. Consequently, its MFNE in Figure 6b is much higher than that of the other algorithms. These boxplots summarize key statistical measures such as quartiles and medians. A smaller inter-quartile range and more compact box shape generally indicate more consistent performance (i.e., higher robustness) across varying elevations. The results suggest that the proposed algorithm exhibits superior consistency compared to the baseline methods in this regard.

6.3. Performance Under Noise

In this subsection, we investigate the impact of observational noise on the algorithm’s reconstruction performance. We evaluate the reconstruction performance of the proposed and baseline schemes for $[eqn]$ values ranging from 2 dB to 22 dB. For these experiments, the down-sampling factor $[eqn]$ was fixed at 10, the MER at 80%, and the number of Monte Carlo trials (M) at 50.

Figure 6e,f illustrate the influence of varying $[eqn]$ levels on the reconstruction performance. Across the tested noise levels, the proposed algorithm consistently demonstrates superior reconstruction performance and exhibits greater resilience to noise compared to the other methods.

To quantitatively assess this noise resilience, we define a metric termed “Slope” as the average absolute rate of change of $[eqn]$ with respect to $[eqn]$ :

$$[eqn]$$

where $[eqn]$ and $[eqn]$ represent the n-th smallest $[eqn]$ value and its corresponding reconstruction error, respectively. $[eqn]$ denotes the total number of distinct $[eqn]$ levels evaluated in the experiment. A smaller Slope value indicates better noise robustness. Table 2 presents the Slope values for all algorithms in both scenarios. The percentage improvement in this Slope metric for the proposed algorithm, relative to the baseline methods, ranges from 7.11% to 95.92%, indicating superior noise adaptability. In scenario $[eqn]$ , the WGDOP method’s low Slope value is notable, primarily because NLOS errors, rather than observational noise, dictate its reconstruction accuracy. Figure 6d shows high reconstruction errors for this method at elevations from 0.5 km to 1.4 km due to multipath effects. Thus, despite good noise robustness, its overall reconstruction performance is substantially inferior to the proposed algorithm, corroborating the conclusions in Table 2.

6.4. Performance with Available Directional Measurements

In this subsection, we evaluate the proposed algorithm’s performance in a scenario where directional measurements of positioning errors are available. This validates the extended algorithm introduced in Section 5.2, assessing its reconstruction accuracy for both scalar ( $[eqn]$ ) and vector ( $[eqn]$ ) representations of the spatial distribution of positioning error.

Table 3 presents performance metrics under varying degrees of MER. Here, $[eqn]$ quantifies the reconstruction error for the $[eqn]$ . $[eqn]$ is computed as the mean of the individual MFNEs obtained for each directional component tensor: $[eqn]$ . Comparing the results in Table 3 with those presented for the scalar error model in Figure 6a, it is found that incorporating directional information enhances the proposed algorithm’s performance. This improvement is observed in the reconstruction of both scalar and vector distributions of positioning errors across various $[eqn]$ conditions.

To more intuitively demonstrate the reconstruction performance for $[eqn]$ , Figure 7 depicts heatmaps of the true and reconstructed directional component tensors in $[eqn]$ and $[eqn]$ . For enhanced visualization, particularly when positioning errors are small, the heatmap values are presented on a dB scale, specifically showing $[eqn]$ for $[eqn]$ .

6.5. Real-Data Experiment

In this subsection, we further validate the proposed algorithm using a real-world dataset collected from a TDOA positioning system. The dataset covers an area of $[eqn]$ m^2^ and includes measurements across 120 frequency bands, ranging from $[eqn]$   $[eqn]$ $[eqn]$ to $[eqn]$   $[eqn]$ $[eqn]$ . The experimental setup comprises four sensors and an unmanned aerial vehicle (UAV) acting as the emitter. The sensors are deployed in a Y-shaped configuration, with baseline distances of 1 $[eqn]$ $[eqn]$ , $[eqn]$   $[eqn]$ $[eqn]$ , and $[eqn]$   $[eqn]$ $[eqn]$ from the peripheral sensors to the central reference sensor, respectively. Figure 8 illustrates the experimental site, including the sensor layout and four sets of flight trajectories for the UAV.

A key challenge presented by this dataset is its extreme sparsity. For example, the trajectory $[eqn]$ is sampled at 386 distinct locations, with ground-truth positions from its onboard GPS. To analyze the spatial distribution of positioning errors, the target area was discretized into a fine-grained $[eqn]$ grid. Consequently, the 386 measurement points represent only a tiny fraction of the 400,000 grid points, resulting in a MER of approximately 99.9%. This high-sparsity scenario serves as a stringent test for the algorithm’s tensor completion capabilities. To evaluate the algorithm’s performance under these conditions, we adopt a cross-validation approach. For each Monte Carlo trial, a subset of the observed data is randomly selected for training the algorithms, with the remainder reserved for testing. The proportion of data used for training is denoted by $[eqn]$ . A total of $[eqn]$ Monte Carlo trials are conducted for each value of $[eqn]$ .

Table 4 presents the percentage of performance improvement of the proposed algorithm compared with the comparative methods in various trajectories. The final column indicates whether an algorithm was capable of reconstructing the complete tensor (“Y” for yes, “N” for no), as some methods could only reconstruct partial slices or fail under extreme sparsity. The performance of the proposed algorithm is comparable to that of the classical GDOP method under simple channel conditions (such as $[eqn]$ ), while it has obvious advantages in complex channel environments (such as $[eqn]$ ). It can also be observed from Table 4 that under relatively ideal environmental conditions ( $[eqn]$ ), the proposed algorithm yields only marginal performance improvements compared to simpler baseline algorithms (e.g., WGDOP), despite possessing significantly higher computational complexity. This raises the question of determining the appropriate scenarios for its application. In practice, however, error map reconstruction is typically performed offline to evaluate the station deployment of positioning systems; thus, time complexity is not the primary concern. In such cases, it is recommended to employ the proposed algorithm across all scenarios, trading higher computational complexity for a performance gain of approximately 30%. Conversely, for time-sensitive scenarios where computational efficiency is critical, we recommend restricting the use of the proposed algorithm to complex terrains, such as urban areas or canyons.

Then, we take trajectory $[eqn]$ as an example to give a further analysis. Table 5 presents the average $[eqn]$ of the reconstructed positioning error tensor on the test sets in trajectory $[eqn]$ , along with computation times.

When comparing the proposed tensor-based modeling approach with WGDOP, we observe that our method further reduces the $[eqn]$ by over 27.96%. This highlights the benefit of our model in capturing complex error distributions. In contrast, traditional data-driven interpolation methods such as RBF and Kriging struggle significantly under such highly sparse conditions, particularly when entire rows or columns (or more generally, large contiguous regions) of the tensor lack measurements, leading to their failure in reconstructing the complete tensor. The NNM-T algorithm was also tested but is omitted from the table as it failed to converge reliably under this level of data scarcity. It also indicates that the proposed algorithm incurs a longer computation time compared to some baselines. However, real-time inference is often not a primary requirement because the analysis of the positioning error map is typically performed offline as a post-processing step.

Figure 9 presents heatmaps of the reconstructed positioning error map from the real-world data using the proposed algorithm. It suggests that the proposed method effectively models underlying error characteristics and can provide valuable guidance for system analysis. For instance, examining the magnitudes of the reconstructed directional error components in Figure 9b–d, we observe that the error magnitudes typically follow the relationship $[eqn]$ . This implies that the positioning reliability is highest along the z-axis and lowest along the x-axis for this specific setup, an observation consistent with statistical analysis of the raw positioning data.

7. Conclusions

This paper tackles the critical challenge of high-fidelity reconstruction of spatial localization error maps from limited measurements. We propose a novel tensor-based framework that overcomes the limitations of conventional model-driven approaches, which are constrained by idealized assumptions, and purely data-driven methods that are dependent on dense observational data. Central to our approach is the TGDOP model, a novel representation that captures the complexity and anisotropic nature of localization errors. To address the practical issue of data scarcity, we developed a physics-informed tensor completion algorithm. By incorporating prior knowledge derived from the analytical properties of the error covariance matrix directly into the factorization process, our algorithm achieves robust reconstruction even from severely incomplete data. Comprehensive experiments on both simulated and real-world TDOA datasets demonstrate the superiority of the proposed method. Specifically, under extremely sparse conditions, our method improves reconstruction accuracy by at least 27.96% compared to state-of-the-art baselines. Future work will explore more efficient optimization strategies for large-scale scenarios and investigate advanced constraints on factor matrices to capture richer physical information.