A Ray-Tracing Based Single-Site Localization Method for Non-Line-of-Sight Environment

1. Introduction

Traditional localization methods assume a LOS condition between the user and the localization station [1]. However, in complex urban environments, there are often intricate NLOS conditions [2]. Obstacles in the LOS path cause reflections, diffractions [3], and transmissions [4] in signal propagation. In such NLOS scenarios, observed time of arrival(TOA) [5], time difference of arrival(TDOA) [6], AOA [7], and received signal strength (RSS) [8] can introduce significant errors, leading to increased inaccuracies in traditional localization methods [9].

Current research on NLOS localization methods primarily focuses on two aspects: NLOS identification techniques and NLOS error mitigation techniques [10,11,12]. NLOS error identification primarily uses empirical measurement data as support, often applied in multi-sensor cooperative localization[13,14,15,16]. Channel information such as RSS, root mean square delay, and root mean square angle spread is used to identify LOS or NLOS nodes[17,18,19,20,21,22]. Some studies also mention using deep learning tools as a key to identifying localization errors. For instance, [23] employed ultra-wideband devices (UWB) and integrated deep learning to develop a deep Q-Learning energy-optimized LOS/NLOS framework, effectively improving the localization accuracy of TDOA algorithms in multi-station scenarios. In [24] , the author utilized an AdaBoost network classifier to identify NLOS nodes, enhancing the TOA localization accuracy. Additionally, Kalman filtering methods are used in some studies to filter out NLOS information [25,26,27]. However, in severe NLOS conditions, the LOS path may not be obtainable, leading to degraded localization accuracy for methods relying on LOS path identification [28]. The second approach is NLOS error mitigation techniques. Existing error mitigation methods mainly rely on data bias or noise deviations in known NLOS environments, which limits the applicability of the algorithms in different NLOS scenarios. In [29], the author assumes all paths are NLOS and introduces a TOA-RSS combined NLOS bias parameter. This method transforms the original non-convex problem into a generalized trust region sub-problem (GTRS) framework, effectively improving localization in NLOS scenarios. In[30], the author employed an equality constrained Taylor series robust least squares method to suppress NLOS errors in UWB localization systems. This approach effectively reduces the algorithm's dependence on NLOS prior data. The work in [31] employed a nonlinear weighted least squares method to eliminate NLOS errors in TOA and RSS hybrid localization. By using a convex hull constraint, this method does not require prior NLOS statistical information. The work in [32] proposed a regulation term least square-based semidefinite programming (RTLS-SDP) method that improves TOA localization accuracy without prior NLOS information. However, these studies address mixed LOS/NLOS scenarios, and the localization performance of these algorithms degrades in purely NLOS conditions.

Additionally, some researchers are focused on using NLOS information for localization. Chee Kiat Seow proposed a method that utilizes single-scattering paths to fuse AOA and TOA for localization [33]. In this approach, multiple scattering paths are identified and removed, which reduces the impact of complex NLOS paths on localization. However, this method does not address localization challenges in severe NLOS environments. In [34], the author first introduced the concept of virtual stations (VSs) and developed the T two-step weighted least squares (TSWLS) method by incorporating TOA observations. This approach achieved localization in simple geometric NLOS scenarios. The work in [35] and [36] proposed a VS localization method based on geometric maps. This method converts NLOS paths into LOS paths and establishes a joint TSWLS equation for TOA and AOA, effectively reducing the impact of some NLOS paths on localization. However, clock synchronization errors among multiple stations can affect localization accuracy. In [37], the author introduces a system of nonlinear equations by establishing a VS and combining AOA and TDOA estimation parameters. The TSWLS method is applied to achieve single-station localization, effectively mitigating the errors caused by clock synchronization in TDOA localization. However, this approach considers only basic geometric information and overlooks the characteristics of electromagnetic wave propagation channels. Consequently, its localization performance is significantly limited in complex NLOS environments.

To address the limitations of traditional localization algorithms in NLOS scenarios, this paper integrates RT algorithms [38,39] with conventional AOA localization algorithms to develop the ray-tracing localization-based service (RT-LBS) method. This approach effectively enables the localization of non-cooperative target source (NCTS) in urban environments using only a single station. Simulations and experimental results demonstrate that the proposed algorithm achieves the CRLB.

In summary, this paper makes the following contributions:

• To improve the performance of existing RT methods and better align them with the proposed localization approach, this paper integrates the advantages of both RT and the shooting and bouncing rays (SBR) technique. Furthermore, an innovative adaptive ray tube structure is introduced, allowing the propagation effects within the environment to be more accurately captured and reflected in the localization algorithm.
• The localization method integrates AOA localization with RT to construct nonlinear equations using GSs generated by sensors in the environment. A heuristic approach for determining equation weights based on angle and power residuals is constructed, and the IRLS method is applied for precise localization of NCTS. Simulations and experimental results show that the proposed method reaches the CRLB.
• In NLOS scenarios, regions with rapid multipath birth and death processes can severely compromise the robustness of localization algorithms. To address this, this paper introduces a MSDM, designed based on multipath similarity. By integrating MSDM, the robustness of localization algorithm in these challenging regions is significantly improved.
• In this paper, the localization algorithm requires repeated invocations of the ray tracing (RT) process for path generation, which involves extensive traversal of the node tree structure constructed by the RT algorithm. This results in decreased computational efficiency. To address this issue, a fast GPU-based algorithm is proposed, which accelerates the path generation component of the RT within the localization algorithm, thereby significantly improving its overall efficiency.

The rest of the paper is organized as follows: Section 2 introduces improvements to the ray tracing algorithm. Section 3 illustrates the construction process of the RT-LBS method. Section 4 showcases the relevant simulations and experimental results supporting the fundamental theory of RT-LBS. Section 5 discusses the findings and limitations of this study. Section 6 provides conclusions and future research plans.

2. The Improvement of the Ray-Tracing Algorithm

RT methods can be divided into two main categories: the image-based ray tracing method (IM-RT) [40] and the shooting and bouncing rays method (SBR) [41]. The IM-RT algorithm first constructs a mirror source tree structure within the environment, followed by solving the path to the receiving point. In contrast, the SBR method emits rays uniformly in all directions from the transmitting point and calculates energy paths by employing a receiving sphere. Due to the differing mechanisms of path generation, the accuracy of electromagnetic field calculations using IM-RT is generally higher than that of SBR. The 2.5D RT method proposed in this paper combines the advantages of IM-RT and SBR. In Figure 1, a simplified flow diagram of the implementation of the RT model is presented. The corresponding process can be divided into the following main steps: simulation setup, loading and preprocessing of environmental information, creation of propagation paths between the Tx and Rx, calculation of fields and the process of simulation results, and the output of simulation results.

The following section provides a detailed explanation of the key components of the proposed algorithm, including the ray splitting algorithm, ray reception algorithm, and electromagnetic field computation. First, the overall algorithm adopts the ray tube launching method, where ray tubes are uniformly launched in a two-dimensional plane, and each ray tube is then traced as it undergoes reflection, transmission, diffraction, and other propagation characteristics in the environment. To further improve accuracy, the proposed algorithm employs ray tube reception instead of the receiving sphere as the reception determination criterion. Finally, the coarse path to the receiving point is refined using the IM method to obtain an accurate path. As shown in Figure 2, the propagation of rays results in a binary tree structure of ray nodes [42].

2.1. Ray Splitting Algorithm

The 2.5D RT algorithm proposed in this paper combines the advantages of the ray tube method and the image method. In the algorithm, rays are designed in the form of ray tubes with splitting attributes. During the initial path search, the uniform splitting threshold of the ray tubes is set based on the size of the geometric facets in the environment (in 2D, this corresponds to the length of line segments). When rays intersect with the environment, the algorithm computes the number of ray splits and generates the corresponding splitting data accordingly. The GS nodes of the split rays are then divided into multiple split GS nodes according to the splitting information. Finally, the algorithm places the split GS nodes into a simulation stack for iterative calculations until all rays are fully traced. Figure 3 presents a schematic diagram of the ray-splitting algorithm, while Algorithm 1 provides the corresponding steps in pseudocode. The coordinates of a GS node $p_{GS}$ are determined by backtracking the propagation distance $l$ from the intersection point $p_{it}$ (e.g., reflection, transmission, diffraction, etc.) between the ray and the environment, in the opposite direction of the outgoing ray $\overrightarrow{r_o}$.

$$p_{GS}=p_{it}-l\cdot \overrightarrow{r_o}.$$

(1)

Algorithm 1. Ray splitting algorithm

Input:              Ray splitting threshold  $t$ Precondition: Create a stack object stack and stack the virtual root node.

While stack is not empty       Pop up the top element node ns on the stack       Obtain the ray r in ns.       Calculate the intersection point between r and scene       Calculate the ray splitting information $\mathsf{\xi }$       $\mathbf{if}\mathsf{\xi }>\mathit{t}$

Get the GS node gns of ns             Calculate the set of splitting nodes $\left\{{\mathit{g}}_{\mathit{i}}\right\}$ based on the $\mathsf{\xi }$             Mount all nodes in $\left\{{\mathit{g}}_{\mathit{i}}\right\}$ onto the right sibling node of ns             delete ns       End if       Trace rays and execute other programs

End While

2.2. Ray Reception Algorithm

In SBR method, the determination of ray reception is based on the concept of a receiving sphere. Due to the uneven distribution of rays in space, it is necessary to dynamically adjust the radius of the receiving sphere to achieve appropriate energy reception. The adaptive splitting ray tube in this algorithm effectively avoids the limitations of the receiving sphere mechanism. As shown in Figure 4, when determining whether point p is received by ray tube$r$, a GS$O^{\prime }$is constructed at point R relative to$O$. If the condition in Equation (2) is satisfied, it can be concluded that point$p$ is captured by ray tube $r$.

$$\left\{\begin{array}{l}{t}_1\le ‖\overrightarrow{O^{\prime }P}‖\le t_1+t_2\hfill \\ \arccos \left({\displaystyle \frac{\overrightarrow{O^{\prime }P}\cdot \overrightarrow{O^{\prime }E}}{‖\overrightarrow{O^{\prime }P}\cdot \overrightarrow{O^{\prime }E}‖}}\right)\le \theta \hfill \end{array}\right..$$

(2)

2.3. Electromagnetic Field Compution

Once the binary tree structure of the ray nodes is established, a depth-first search (DFS) traversal is conducted on the ray tree, yielding multiple potential paths from the transmitter (Tx) to the receiver (Rx). If there is a LOS condition between the Tx and Rx, the radiation field at the arrival point p is calculated according to Equation (3):

$${\stackrel{\rightharpoonup }{E}}_{Los}={\stackrel{\rightharpoonup }{E}}_0{\displaystyle \frac{e^{-jk{r}_0}}{r_0}},$$

(3)

where${\stackrel{\rightharpoonup }{E}}_0$denotes the electric field strength at the position of the Tx antenna, k stands for the wave number, and$r_0$ is the total length of LOS propagation path.

In the case where a NLOS exists between the Tx and Rx, the LOS field ${\stackrel{\rightharpoonup }{E}}_0$ is first computed using Equation (4), and then the contributions of other electromagnetic propagation mechanisms to the field strength, ${\stackrel{\rightharpoonup }{E}}_{inc}$, is progressively calculated along the ray path.

$${\stackrel{\rightharpoonup }{E}}_{NLOS}={\stackrel{\rightharpoonup }{E}}_{inc}\cdot {\displaystyle \prod }_{h=1}^n\overline{\overline{R_h}}\cdot {\displaystyle \prod }_{i=1}^m\overline{\overline{D_i}}\cdot {\displaystyle \prod }_{t=1}^l\overline{\overline{T_t}}\cdot {\displaystyle \prod }_{s=1}^{n+m+l}{A}_s\cdot {\displaystyle \prod }_{q=1}^{n+m+l}{e}^{-jk{r}_q},$$

(4)

where n represents the number of reflections, m represents the number of diffractions, and t represents the number of transmissions along the path. $A_s$ is the spherical wave spreading factor, and $r_q$ is the distance between the q-th node and the (q+1)-th node.

The total field strength at the receiving point is obtained by combining the field contributions from all paths, which is the sum of the field strengths of all rays. It should be noted that this paper presents two modes for combining the field strengths: $E_{total}$, which is the scalar combination of fields, and ${\stackrel{\rightharpoonup }{E}}_{total}$, which is the vector combination of fields.

$$\left\{\begin{array}{c}{E}_{total}={\displaystyle \sum _{i=1}^l\left|\stackrel{\rightharpoonup }{E_i}\right|}\\ {\stackrel{\rightharpoonup }{E}}_{total}={\displaystyle \sum _{i=1}^l\stackrel{\rightharpoonup }{E_i}}\end{array}\right.,$$

(5)

where $l$ denotes the number of the ray paths.

After deriving the field strength formula, by incorporating the Tx antenna gain $G_t$, the receiver antenna gain $G_r$, and the antenna radiated power $P_t$, the received power $P_r$ at the field point can be obtained as follows:

$$\begin{array}{c}{P}_r=10\lg \left(P_{tx}{G}_{tx}{G}_{rx}{\left({\displaystyle \frac{\lambda }{4\pi }}\right)}^2{\left|{\displaystyle \frac{E_{total}}{E_0}}\right|}^2\right)+30 dBm\end{array},$$

(6)

3. Single-Site Localization Algorithm for NLOS Environment

The RT-LBS method proposed in this paper is designed for localization in complex geometrically obstructed environments. This work effectively utilizes the RT method's in section 2 predictive capabilities in complex scenarios, extending traditional AOA localization algorithms. The proposed algorithm enables accurate NLOS localization with a single station. A simplified flowchart of the RT-LBS algorithm is provided in Figure 5. The associated process can be divided into the following key steps: establishment of generalized sources, generalized sources filtering rules, generalized sources weight calculation, and linear equation solving.

3.1. Establishment of Generalized Sources

The RT algorithm can track each ray emitted from the source, thus covering most of the geometric surface features within the environment. In the localization phase, the sensor is treated as a virtual source. This sensor is assumed to receive a set of AOAs $\left\{{\varphi }_1,{\varphi }_2,\dots ,{\varphi }_n\right\}$ corresponding to the electromagnetic waves emitted by a NCTS, where n represents the number of angle information received by the sensor. Furthermore, it is assumed that the angle measurement error of the sensor follows a zero-mean Gaussian distribution, with the standard deviation${\sigma }_{\varphi }$, Therefore, each angle interval received by the sensor is $\left[{\varphi }_i-{\sigma }_{\varphi },{\varphi }_i+{\sigma }_{\varphi }\right]$. Next, ray tubes are uniformly emitted within the angle interval at increments of $\Delta \varphi$. The algorithm tracks each ray tube and constructs a binary node tree structure. By traversing each node of the binary tree structure, all GSs related to the sensor are obtained.

The position of the reflected GS $GSr$ is the location of the mirror image of source $S$ in the geometric environment. The position of the diffracted GS $GSd$ is located on the diffraction wedge that generates the diffracted rays. The position of the transmitted GS $GSt$ is obtained by backtracking the propagation distance $t$ from the transmission node location along the outgoing transmission direction. It is important to note that the involved reflected, transmitted, and diffracted GSs are not limited to first-order. They can undergo multiple reflections, transmissions, and diffractions, forming higher-order GSs through combinations of these interactions. By traversing all possible propagation path nodes in the environment, all GSs related to $S$ are obtained. This allows the corresponding AOA localization equations to be formulated as follows:

$$\left\{\begin{array}{c}\left(x-x_1\right)\tan {\varphi }_1=y-y_1\\ \left(x-x_2\right)\tan {\varphi }_2=y-y_2\\ ⋮\\ \left(x-x_m\right)\tan {\varphi }_m=y-y_m\end{array}\right.,$$

(7)

where m denotes the number of the GS, $\varphi$ stands for the AOA received by GS.

3.2. Generalized Sources Filtering Rules

The RT algorithm generates a large number of GSs, each carrying angle information about the NCTS, which can form m localization equations based on AOA. There is exactly one valid GS in each set of GSs for a given AOA that is in a LOS environment with the NCTS (here, we assume that there is only one multipath arriving at the receiving sensor per angle).

Therefore, the system of linear equations formed by the GS set contains a large amount of noise, making the system ill-conditioned, which in turn causes the target localization solution to fail to converge. To address this issue, this paper proposes a filtering method specifically designed for the GS set. This method effectively eliminates invalid equations, improving the stability and accuracy of the localization solution.

In the first step, by selecting two GSs from the GS set and combining them $Pai{r}_{i,j}=\left\{G{S}_i,G{S}_j\right\},i\ne j$, a system of equations can be constructed based on the two selected GSs:

$$\left\{\begin{array}{c}\left(x-x_i\right)\tan {\varphi }_i=y-y_i\\ \left(x-x_j\right)\tan {\varphi }_j=y-y_j\end{array}\right..$$

(8)

This equation has an analytical solution, denote as $P_{i,j}=\left(x_p,y_p\right)$:

$$\left\{\begin{array}{c}{x}_p={\displaystyle \frac{\left(x_i\sin {\varphi }_i-y_i\cos {\varphi }_i\right)\cos {\varphi }_j-\left(x_j\sin {\varphi }_j-y_j\cos {\varphi }_j\right)\cos {\varphi }_i}{\sin {\varphi }_j}}\\ y_p={\displaystyle \frac{\left(x_i\sin {\varphi }_i-y_i\cos {\varphi }_i\right)\sin {\varphi }_j-\left(x_j\sin {\varphi }_j-y_j\cos {\varphi }_j\right)\sin {\varphi }_i}{\sin {\varphi }_j}}\end{array}\right..$$

(9)

The GSs filtering rules proposed in this paper are as follows:

• If $P_{i,j}$ is located inside a building or outside the solution domain, the position is considered invalid;
• Construct line segments $s_1$ and $s_2$ connecting position $P_{i,j}$ with $G{S}_i$ and $G{S}_j$ at the intersection points $nod{e}_i$ and $nod{e}_j$, respectively. If $s_1$ and $s_2$ intersect within the environment, then $P_{i,j}$ is considered an invalid solution.

In the second step, begin looping through the $C_m^2$ pairs of GSs, solving for the validity of each GS pair. If $Pai{r}_{i,j}$ is valid, increment the count weight $w$ for the corresponding GSs by 1. After traversing all GSs, delete those with a count weight $w=0$.

Additionally, this paper employs a hash mapping method to eliminate duplicate GSs with overlapping positions. During the hash mapping process, angle data for the identified duplicate GSs are averaged. However, GSs generated by diffraction and transmission nodes are not subject to this duplication removal.

3.3. Generalized Sources Weight Calculation

Let the remaining number of GSs be $l$. By combining these remaining sources, $C_l^2$ pairs of GS combinations are generated, with the solution for each pair denoted as$p_{i,j}$. The solution set $\left\{p_{i,j}\right\}$ from the $C_l^2$ combinations is then clustered based on a given inter-cluster distance threshold $d$, resulting in $q$ clusters. The center position of each cluster is denoted as $p_i,i\in \left[1,q\right]$.

In this study, a generalized residual calculation method is proposed. The method first calculates the path from $p_{i,j}$ to the sensor $S$, and then, using electromagnetic computation techniques, determines the angles and power of the $n_1$ multipath components arriving at $S$. These angles and power values are subsequently sorted in descending order based on power magnitude. Since the power of the NCTS is unknown, it is not feasible to directly compute the residuals based on absolute power values. To address this problem, a power difference of arrival (PDOA)-based residual calculation method is proposed in this study. The original RSS at the sensor is denoted as $p{r}_0$, while the calculated RSS at the sensor from the NCTS​ is denoted as$p{r}_1$. The PDOA is then computed, and combined with Equations (4) and (6), we can obtain:

$$\begin{array}{l}\Delta pr=p{r}_1-p{r}_0\\ =10\lg \left(P_{tx}{G}_{tx}{G}_{rx}{\left({\displaystyle \frac{\lambda }{4\pi }}\right)}^2{\left|{\displaystyle \frac{E_{total}^{\left(1\right)}}{E_0}}\right|}^2\right)+30-\left(10\lg \left(P_{tx}{G}_{tx}{G}_{rx}{\left({\displaystyle \frac{\lambda }{4\pi }}\right)}^2{\left|{\displaystyle \frac{E_{total}^{\left(0\right)}}{E_0}}\right|}^2\right)+30\right)\\ =20\lg \left({\displaystyle \frac{E_{total}^{\left(1\right)}}{E_{total}^{\left(0\right)}}}\right)\\ =20\lg \left({\displaystyle \frac{{\displaystyle \sum {\displaystyle \prod }_{h=1}^{n_1}\overline{\overline{R_h}}\cdot {\displaystyle \prod }_{i=1}^{m_1}\overline{\overline{D_i}}\cdot {\displaystyle \prod }_{t=1}^{l_1}\overline{\overline{T_t}}\cdot {\displaystyle \prod }_{s=1}^{n_1+m_1+l_1}{A}_s\cdot {\displaystyle \prod }_{q=1}^{n_1+m_1+l_1}{e}^{jk{r}_q}}}{{\displaystyle \sum {\displaystyle \prod }_{h=1}^{n_0}\overline{\overline{R_h}}\cdot {\displaystyle \prod }_{i=1}^{m_0}\overline{\overline{D_i}}\cdot {\displaystyle \prod }_{t=1}^{l_0}\overline{\overline{T_t}}\cdot {\displaystyle \prod }_{s=1}^{n_0+m_0+l_0}{A}_s\cdot {\displaystyle \prod }_{q=1}^{n_0+m_0+l_0}{e}^{jk{r}_q}}}}\right)\end{array}.$$

(10)

The definitions of the variables in the equation can be found in Equations (4) and Equation (6) . From Equation (10), it can be derived that $\Delta pr$ is solely related to the number and coefficients of reflections, diffractions, and transmissions between the two paths, and is independent of the initial field. Therefore, $\Delta pr$, can be effectively used to quantify the residual between two positions.

Therefore, after obtaining the original angle-power data $\left\{{\varphi }_{raw,i},\Delta p{r}_{raw,i}\right\},i\in \left[1,n\right]$ and the simulated angle-power data $\left\{{\varphi }_{sim,i},\Delta p{r}_{sim,i}\right\},i\in \left[1,n_1\right]$, it is often the case that the number of multipath components in the measured and simulated data are not identical, i.e., $n_1\ne n$. As a result, it is not possible to directly compute the residuals. This paper proposes a generalized residual computation method, which first employs the Hungarian algorithm (as provided by the GLPK open-source library) to determine the minimum residual $\arg \min \left\{{\varphi }_{raw,i},{\varphi }_{sim,j}\right\}$, between two sets of angle data. This process yields the angular residual, $r_{\varphi }$ and the corresponding assignment sequence. Using the assignment sequence, the power difference residual $r_{\Delta pr}$ is calculated. In cases where $n_1<n$, the method also records the number $k=n-n_1$ of unmatched original data points.

$$\left\{\begin{array}{l}{r}_{\varphi }={\displaystyle \sum _{i=1}^{\min (n,n_1)}{‖{\varphi }_{raw,i}-{\varphi }_{sim,j}‖}_2}\hfill \\ r_{\Delta pr}={\displaystyle \sum _{i=1}^{\min (n,n_1)}{‖\Delta p{r}_{raw,i}-\Delta p{r}_{sim,j}‖}_2}\hfill \end{array}\right.,$$

(11)

After calculating the angular residual $r_{\varphi ,i}$ and the power difference residual $r_{\Delta p,i}$ for all clusters, the maximum values $r_{\varphi ,\max }$ and $r_{\Delta p,\max }$, as well as the average value ${\overline{r}}_{\varphi }$ and ${\overline{r}}_{\Delta p}$ of the residuals are determined. Given that instances of unmatched points arise during the residual calculation process, the total residual needs to be adjusted by adding k-times the average residual, as expressed by the following equation:

$$\left\{\begin{array}{c}{r}_{\varphi ,i}=r_{\varphi ,i}+k{\overline{r}}_{\varphi }\\ r_{\Delta p,i}=r_{\Delta p,i}+k{\overline{r}}_{\Delta p}\end{array}\right..$$

(12)

Subsequently, the residuals for cluster are normalized based on the maximum values of residuals $r_{\varphi ,\max }$ and $r_{\Delta p,\max }$, yielding the normalized residuals ${\widehat{r}}_{\varphi ,i}$ and ${\widehat{r}}_{\Delta p,i}$. In addition, the size of the cluster is incorporated into the weight calculation, with the normalized weight of the cluster denoted as ${\widehat{w}}_{cluster,i}$. The weight expression is then derived as shown:

$$w_i=a{\displaystyle \frac{1}{{\widehat{r}}_{\varphi ,i}+\epsilon }}+b{\displaystyle \frac{1}{{\widehat{r}}_{\Delta p,i}+\epsilon }}+c{\widehat{w}}_{cluster,i},$$

(13)

where a, b, and c represent the coefficients of the AOA weight, PDOA difference weight, and cluster size weight, respectively. These coefficients must satisfy condition $a+b+c=1$, which is typically set to the same weight. $\epsilon$ is a small decimal value, typically set to 1e−6.

After determining the weight of each cluster, the weights of the individual GSs within the clusters are updated, yielding a system of linear equations with the assigned initial weights.

3.4. IRLS

After determining the weights for each GS, an IRLS method based on the initial weights is proposed to solve the system of linear equations. The AOA localization equation can be reformulated in terms of residuals as follows:

$$r_i\left(x,y\right)=\left(x-x_i\right)\tan {\varphi }_i-\left(y-y_i\right).$$

(14)

The objective function $S\left(x,y\right)$ to be optimized in the IRLS method is the weighted sum of squared residuals, which is the squared ℓ2​-norm of the residuals, weighted as follows:

$$S\left(x,y\right)=W{‖r‖}_2^2={\displaystyle \sum _{i=1}^m{w}_i{\left(r_i\left(x,y\right)\right)}^2},$$

(15)

where $w_i$ denotes the weight of each equation.

In the IRLS iteration process, the initial position $\left(x_0,y_0\right)$ of the NCTS is first estimated. During each subsequent iteration, the residual $r_i\left(x_i,y_i\right)$ for each equation is calculated, and the weights are updated accordingly:

$$w_{i+1}={\displaystyle \frac{1}{{\left(r_i\left(x_i,y_i\right)\right)}^2+\epsilon }},$$

(16)

The estimated position after the k-th iteration can be obtained based on the least squares method as follows:

$$\left[\begin{array}{c}{x}_{i+1}\\ y_{i+1}\end{array}\right]=\left[\begin{array}{c}{x}_i\\ y_i\end{array}\right]-{\left(J^{T}WJ\right)}^{-1}{J}^{T}Wr,$$

(17)

where $J$ represents the Jacobian matrix of the residual equation, and $r$ denotes the residual vector.

Additionally, an error threshold must be established. Typically, the iteration stopping criterion can be defined as:

• Whether the value of the objective function $S\left(x,y\right)$ falls below the predefined threshold $s_a$;
• Whether the Euclidean distance between consecutive iterative solutions is less than the threshold $l_a$;
• Whether the number of iterations has reached the maximum limit $k$.

3.5. Cramér-Rao Lower Bound

The CRLB is widely used to evaluate the theoretical limit of localization accuracy under the influence of unbiased noise. It defines the lower bound of localization error in the presence of noise. For the AOA localization method, it is assumed that the measurement angle errors follow a zero-mean Gaussian distribution.

$${\varphi }_i={\tan }^{-1}\left({\displaystyle \frac{y-x_i}{x-y_i}}\right)+{\eta }_i,$$

(18)

where ${\eta }_i$ represents the error term.

The Jacobian matrix is obtained by taking the partial derivatives with respect to x and y, respectively:

$$J_{AOA}=\left[\begin{array}{cc}{\displaystyle \frac{\partial {\varphi }_1}{\partial x}}& {\displaystyle \frac{\partial {\varphi }_1}{\partial y}}\\ {\displaystyle \frac{\partial {\varphi }_2}{\partial x}}& {\displaystyle \frac{\partial {\varphi }_2}{\partial y}}\\ ⋮& ⋮\\ {\displaystyle \frac{\partial {\varphi }_m}{\partial x}}& {\displaystyle \frac{\partial {\varphi }_m}{\partial y}}\end{array}\right].$$

(19)

In accordance with the definition of the Fisher Information Matrix (FIM), the following expression can be derived:

$$I\left(\mathrm{x}\right)={\displaystyle \frac{1}{{\sigma }^2}}{J}^{T}J.$$

(20)

Thus, the CRLB is given by the inverse of the FIM, expressed as:

$$CRLB\left(\mathrm{x}\right)=I{\left(\mathrm{x}\right)}^{-1}.$$

(21)

4. Simulation and Experimental Results

4.1. Measurement Compaign

4.1.1. Measurement Equipment

The architecture of the power measurement system and the key instrumentation are illustrated in Figure 6. A continuous wave (CW) signal is generated by the signal generator and amplified up to a maximum of 43 dBm by the signal amplifier. On the Rx platform, the signal is directly fed from the receiving antenna into a spectrum analyzer, which is controlled by a computer to capture power information. Additionally, a real-time kinematic (RTK) system and its antenna are mounted on a cart to collect real-time positional data. Both the transmitting and receiving antennas are vertically polarized omnidirectional antennas, with a height of 1.85 m for each.

The architecture of the localization test system and the key instruments are shown in Figure 7. The NCTS platform uses the same architecture as the transmission platform in the power measurement system. An 8-channel uniform circular array (UCA) antenna platform is employed to measure the AOA of incoming signals, and the multiple signal classification (MUSIC) algorithm is used to estimate the AOA. Figure 7 presents an image of the 8-channel circular array antenna along with the schematic of the corresponding RF board. Four ADRV9008 chips are used to form the 8-channel receiving system, with clock synchronization achieved through the AD9528 chip. Under anechoic chamber measurement conditions, the direction-finding accuracy can reach approximately 1°.

4.1.2. Measurement Scenarios

The measurements were conducted in Xi'an, Shaanxi Province, China. Within a modern technology park comprising 23 buildings. The area spans 0.2 square kilometers. A geometric model of the technology park was generated using a handheld LiDAR device, which offers an accuracy of up to 5 cm. Figure 8 presents the raw laser point cloud data alongside the reconstructed site floor plan. All measurement scenarios are situated within an urban environment, with the scenarios in this study divided into fading measurement scenarios and localization test scenarios.

During the mobile power measurement process, it is essential to maintain a constant receiver speed of 5 km/h. The measurement duration for each frequency point ranges from approximately 450 to 500 seconds, with the longest measurement path extending up to 700 meters. The street canyons in the measurement area are notably narrow, varying between 6 and 12 meters in width. Throughout the measurement campaign, most of the receiver's path is obstructed by buildings, resulting in NLOS conditions. As the path progresses, the prevalence of these NLOS conditions increases, leading to a substantial rise in propagation losses.

4.2. Processing Measured Results and Verification of Models

4.2.1. Power Measurement Data Processing

A total of five frequency points were measured at 3 GHz, 3.6 GHz, 4 GHz, and 5 GHz during the power measurement process. Figure 9 illustrates the measurement paths for each frequency point, along with the RSS at each location, depicted using a pseudo color map. The Tx was positioned within a street canyon, and most of the receiver's path remained in NLOS conditions relative to the Tx. The detailed measurement parameters are provided in Table 1. It is important to note that a fully automated measurement methodology was adopted, with the RTK system and spectrum analyzer recording real-time positional information (with centimeter-level accuracy) and power at a rate of 20 frames per second. The location and power data were synchronized by matching the timestamps of the positional records with those of the power measurements.

Given the positional inaccuracies in the raw power measurement data, a dB-based sliding window method was applied to process the received data. This approach mitigates power fluctuations caused by phase variations due to minor positional errors. The sliding window size was set to 5 data points. Figure 10 presents the results after applying the sliding filter for the five frequency points, where the black line represents the raw measurement data, and the blue line represents the filtered data.

4.2.2. Power Simulation and Verification

Both the algorithm presented in [39] and the algorithm proposed in this paper utilize a combination of three reflections and two diffractions for the calculations. Table 2 lists the material properties used in the algorithms, including the conductivity $\sigma$ and the relative dielectric constant $\mu$. In the measurement environment of this study, the building walls are composed of glass, the exterior wall decorations are made of aluminum, and the ground material is concrete.

The mean and standard deviation of the comparison between the proposed algorithm and that of reference [39] against the measured data are summarized in Table 3. As illustrated by Figure 11 and Table 3, the proposed algorithm provides a notably superior performance in power prediction compared to the algorithm presented in [39]. For the initial 2500 data points, both algorithms exhibit strong agreement with the measured data. This is primarily due to the fact that, in regions with mild shadowing, the primary path energy is largely dominated by low-order reflections and diffractions, resulting in stable path search performance for both algorithms. However, under severe NLOS conditions (e.g., data points 4000–4600 for 3 GHz in Figure 11), the main path energy becomes more scattered, and the algorithm from reference [39] fails to capture certain multipath components, leading to a mismatch with the measured data and a significant loss of accuracy in deep shadowed areas. In contrast, the proposed algorithm is capable of identifying a greater number of multipath components even in such challenging deep shadow conditions, thereby achieving enhanced computational accuracy in NLOS scenarios.

4.2.3. AOA Measurement Data Processing

The angle measurement scenario was selected in a typical NLOS T-shaped street canyon within the park. The sensor was positioned in a canyon with a width of 12 meters, while three source locations were placed in a canyon with a width of 9 meters. T1 and T3 were located at the two ends of the canyon, both maintaining an NLOS condition with respect to R. T2 was placed at the intersection of the canyon, where it had a LOS condition with respect to R.

Figure 12. Angle measurement scenario and the positions of the NCTS (denotes T1, T2, and T3) and sensor (denotes R).

Figure 12. Angle measurement scenario and the positions of the NCTS (denotes T1, T2, and T3) and sensor (denotes R).

Figure 13, Figure 14 and Figure 15 illustrates the results of 500 AOA measurements conducted at the three sites, T1, T2, and T3, under 5GHz frequency. During the measurement process, both the direction-finding sensor and the source were maintained at a height of 1.85 meters and remained stationary relative to each other. Table 4 summarize the multipath angle power differences extracted from the measured AOA spectra at locations T1, T2, and T3.

4.2.4. Localization Algorithm Verification

Figure 13, Figure 14 and Figure 15 depicts the multipath clusters with higher received power along with their corresponding propagation mechanisms. Figure 16 presents a comparison between the average AOA spectrum derived from 500 measurements and the results obtained using the proposed RT algorithm. The red curve represents the AOA spectrum received at point R, the solid blue line indicates the multipath clusters with energy within 40 dB of the first arrival path, while the light blue line represents the clusters with energy lower than 40 dB of the first path. Table 5 provides the differences between the root mean square (RMS) angular spread (AS) predicted by the proposed algorithm and the measured angular spread. Both the figures and tables indicate that the proposed algorithm demonstrates strong consistency with the measured AOA spectrum.

The accuracy of the proposed localization algorithm was evaluated and verified by combining the measured angles and power values listed in Table 4. Based on the proposed algorithm, the positioning errors for T1, T2, and T3 are documented in Table 6. All positioning errors are within 1m range, with an average error of approximately 0.69 meters.

5. Efficiency and Performance Analysis of Localization Algorithm

5.1. Comparison of Localization Accuracy Under Different AOA and PDOA Errors

In this section, a representative NLOS scenario is employed to evaluate the localization accuracy of the proposed algorithm. The scenario consists of six buildings, each approximately 10 meters in height, spanning an area of 110 m × 150 m. The sensor is located at point R with coordinates (45, 21), while 3 NCTSs A, B, and C are positioned at (76, 56), (121, 74), and (70, 90), respectively. Figure 17 depicts the simulation environment along with the locations of points R (receiver) and A, B, and C (transmitters). The 3 NCTSs (A, B, and C) are all in NLOS conditions relative to sensor R, with the degree of NLOS conditions ranked by severity. The multipath propagation between R and the stations A, B, and C is shown in Figure 17, where solid lines represent multipath clusters with received power within 30 dB of the first path, and dashed lines denote clusters with received power more than 30 dB below the first path.

Single-station localization is performed for each of the points A, B, and C. For point A, the received multipath is predominantly composed of paths involving up to 3 reflections. For point B, the received multipath consists mainly of paths involving up to 4 reflections and one diffraction. For point C, the received multipath comprises paths involving 2 to 4 reflections and one diffraction.

In the first round of simulations, the propagation mechanisms involved in localization include both reflection and diffraction. The sources are positioned at points A, B, and C, with errors due to received power being neglected. The AOA error is assumed to follow a zero-mean distribution with standard deviations ranging from 0.1° to 6°. For each error value, 10,000 simulations were performed. Figure 18 presents a comparison of the localization errors at points A, B, and C against the CRLB. The results demonstrate that the accuracy of the proposed RT-LBS algorithm approaches the CRLB. Moreover, considering the number of multipath components (MPCs) and the propagation distance between points A, B, C, and R, it can be deduced that point A exhibits the highest localization accuracy due to the larger number of MPCs and its relatively shorter propagation distance. In contrast, point C shows the lowest localization accuracy, attributable to its fewer MPCs and longer propagation distance.

In the second round of simulations, the localization accuracy of the algorithm is evaluated under the combined influence of both angle and power errors. The RMSE is calculated based on 10,000 independent simulations. As in the first simulation, the AOA error is modeled as a Gaussian random variable with zero mean and a standard deviation ranging from 0.1° to 6°. The RSS deviation, which directly affects the received multipath power, is also modeled as a Gaussian random variable with zero mean and a standard deviation between 0 dB and 10 dB. Figure 19 illustrates the localization accuracy at point A under different AOA and PDOA error conditions. The proposed algorithm demonstrates that, as PDOA error increases, localization accuracy is impacted to some extent. However, under low AOA error conditions, even with PDOA errors reaching up to 10 dB, the localization accuracy remains relatively high. In contrast, under high AOA error conditions, the localization accuracy degrades linearly with increasing PDOA error. Figure 20 and Figure 21 present the localization errors at points B and C, respectively, for varying AOA and PDOA errors.

The simulation results at points A, B, and C demonstrate that the proposed algorithm is capable of accurately localizing the source not only in environments with low-order reflections and diffractions (e.g., point A) but also in more complex scenarios involving higher-order reflections (e.g., points B and C). Given the prediction accuracy of the ray tracing algorithm proposed in this study, which is approximately 6 dB, the localization accuracy at point A can be maintained within 2 meters when the angle error is around 1°. Similarly, at points B and C, the accuracy remains within 5 meters under the same conditions. These findings highlight the robustness of the proposed algorithm across a range of propagation environments and error conditions.

5.2. Localization Error Analysis Across the Entire Plane

In this section, a study of the localization error across the plane is conducted based on the scenario shown in Figure 17. The direction-finding sensor is placed at point R, while the NCTSs are uniformly distributed over a rectangular area (120m × 80m) with a spacing of 0.5 meters. Both the sensor and the NCTSs are equipped with omnidirectional antennas. The power of the non-cooperative sources is set to 1.0 W, and the operating frequency is 5 GHz. A combination of up to 4 reflections and one diffraction is configured for the GSs search.

As the proposed algorithm is fundamentally based on the ray tracing technique, the accuracy of the localization process is directly contingent upon the precision of the ray tracing algorithm. However, during planar localization, as the source position shifts, certain MPCs may emerge or disappear. This dynamic evolution of MPCs gives rise to regions characterized by significant multipath variability, where the multipath features of adjacent coordinate points exhibit substantial disparities. When the source is located within such regions of pronounced multipath fluctuations, the increase in angular error may lead to a corresponding escalation in localization error. This phenomenon, wherein rapid changes in multipath characteristics result in localization inaccuracies, is referred to as multipath fluctuation error (MFE).

Furthermore, multipath similarity constitutes another critical factor that can contribute to localization error, which denotes multipath similarity error (MSE). In this paper, multipath similarity is defined as the condition where the objects encountered at each propagation node across different MPCs are identical, and the AOA and power levels exhibit high degrees of similarity. In regions where multipath similarity prevails, the localization algorithm is likely to yield similar solutions, thereby linking the localization error to the extent of the similarity region.

Thus, in planar localization, both MFE and MSE emerge as significant sources of degradation in the accuracy of the localization algorithm. These factors underscore the challenges posed by complex multipath environments in achieving high-precision localization.

To address the challenges posed by MFEs and MSEs, this paper introduces a displacement correction method based on a multipath similarity displacement matrix (MSDM). Figure 22 (a) - (f) illustrate the MSDM under AOA error deviations ranging from 0.1° to 6°. Each element in these matrices represents the average adjacent multipath similarity distance (MSD) at the corresponding coordinate point. For instance, if the average MSD at point p is denoted as d, it can be inferred that the MPCs within a circular region centered at point p, with a radius of d, exhibit similar characteristics. As the angular deviation increases, the average adjacent MSD in the matrix also expands, indicating that the range of multipath similarity grows progressively larger. This relationship underscores the increasing spatial correlation of MPCs as angular errors rise.

The displacement compensation expansion method expands the single coordinate solution in the GS cluster to 9 coordinate solutions, as illustrated in Figure 23. To account for multipath fluctuation errors, the algorithm compensates by extending the solution by minimum resolution distance around the coordinate point. For multipath similarity errors, the algorithm expands the solution by the average adjacent distance $d$ from the corresponding matrix around the coordinate point, serving as a compensatory solution.

Figure 24, Figure 25, Figure 26, Figure 27, Figure 28 and Figure 29 illustrate the planar localization accuracy distribution for both the original algorithm and the displacement correction algorithm. The optimized localization algorithm significantly improves accuracy, particularly at the shadow boundaries. The blank areas in the figures represent regions with no solution, primarily due to an insufficient number of multipath signals reaching these areas.

Table 7 summarizes the key performance metrics of the localization algorithm with the displacement correction method applied. The average localization accuracy of the algorithm improves by approximately 0.5 meters, with the proportion of localization errors less than 10 meters increasing by an average of 1.4%.

5.3. GPU Accelerated Analysis of Localization Algorithm

The most computationally intensive component of the proposed localization algorithm lies in the path-finding process within the ray tracing algorithm. Consequently, the acceleration method introduced in this work is specifically designed to optimize the path-finding procedure. As illustrated in Figure 30, the technical workflow of the proposed acceleration method proceeds as follows: First, the CPU performs a DFS of the binary tree structure, storing the resulting path node sequences in GPU memory. Next, the path node capture kernel function is employed to compute the valid path node sequences for the captured receiver points. Finally, the path-finding kernel function is invoked to search for the paths corresponding to each valid node and to apply the necessary corrections. This approach significantly enhances computational efficiency while maintaining path accuracy.

To assess the computational efficiency of the proposed algorithm, simulations were performed using a detailed electronic map of Xi’an, Shaanxi Province, China, extracted from OpenStreetMap. This dataset comprises 6,189 buildings and more than 35,700 facets. The transmitter was positioned at the center of the map, elevated to a height of 200 meters, with a transmission power set to 50 W and an operating frequency of 1 GHz. The receivers were placed 2 meters above ground level and uniformly distributed across a rectangular area of 2 km × 2 km, with a spacing of 1 meter between them. Both the transmitter and receivers employed vertically polarized omnidirectional antennas.

The proposed algorithm is fully implemented in C++/CUDA C++. Simulations were conducted on a desktop computer equipped with an Intel(R) Core(TM) i9-13900K CPU, an NVIDIA RTX 3090 GPU, and 128 GB of memory. To evaluate the computational efficiency of the GPU-accelerated algorithm, partial surface simulations were carried out for the scenario depicted in Figure 31, covering a total of 35,688 points. The simulation parameters were configured to allow up to 3 reflections and 1 diffraction. The simulation runtimes are presented in Figure 32 and Table 8, demonstrating that the GPU-based approach achieved a remarkable speedup of approximately 4,835.6 times compared to the single-core CPU implementation, with an average computation time of 7 millisecond per receiving point.

6. Discussion and Future Work

Based on the series of empirical measurements and simulation tests presented in this paper, the effectiveness and robustness of the proposed NLOS localization algorithm have been clearly demonstrated. However, we identify two primary limitations of the proposed algorithm that warrant further exploration in future research:

(1) The geometric modeling in this study is derived from high-precision LiDAR point cloud data, achieving an accuracy of up to 5 cm. However, this level of precision is minimal compared to the errors induced by multipath propagation, and thus, the potential localization bias resulting from inaccuracies in building structures within the electronic map is not reflected in this work. In practical engineering applications, high-precision point cloud-based modeling is rarely available; instead, digital maps are typically sourced from platforms such as OpenStreetMap or other publicly accessible databases. As a result, in future research, it will be essential to account for the errors arising from the lower accuracy of digital map modeling to enhance the algorithm's adaptability and overall localization performance in practical scenarios.

(2) The localization algorithm proposed in this paper is fundamentally based on the ray tracing algorithm, whereby the accuracy of the RT predictions directly influences the algorithm's localization performance. However, one limitation of the current model is that it does not account for the impact of vegetation penetration on the propagation model, which leads to a decrease in the predictive accuracy of the RT algorithm in areas affected by vegetation. To address this issue, future work should incorporate a vegetation penetration model into the current propagation framework to enhance the overall accuracy of the localization predictions in such environments.

Furthermore, we suggest the following areas as key focal points for future research:

(1) In this paper, we proposed a localization algorithm that relies solely on a single direction-finding sensor. However, many modern electronic devices are equipped with more than one type of sensor, such as Wi-Fi, Bluetooth, and UWB, all of which possess direction-finding and ranging capabilities. Therefore, incorporating AOA with other techniques, such as TOA and TDOA, could serve as a hybrid approach to mitigate the limitations associated with high AOA errors.

(2) In this paper, the localization algorithm was evaluated in urban NLOS scenarios. However, NLOS conditions are often more prominent in indoor environments, suggesting that the algorithm's performance in indoor localization warrants further investigation. Moreover, in enclosed environments such as indoors, specific propagation mechanisms involving floors and ceilings limit the effectiveness of the current 2D RT algorithm. This limitation hinders its ability to accurately trace ray paths interacting with planar structures like floors and ceilings in enclosed or semi-enclosed spaces. Thus, developing a 3D localization algorithm holds significant potential for improving accuracy in such scenarios.

7. Conclusions

This paper proposes a single-station localization method that integrates a RT algorithm specifically designed for NLOS environments. In complex NLOS scenarios, signals emitted by the source undergo multipath propagation before reaching the sensor. The proposed method exploits the AOA and the received power differences to construct a GS set with initial weights, based on a GS filtering rule. Each GS corresponds to a linear equation derived from the AOA. By incorporating the IRLS method, the algorithm achieves high-accuracy localization of the target. Simulations and experiments demonstrate that with an AOA error of 1°, the outdoor localization accuracy can reach approximately 1 meter. A key advantage of the proposed approach is that it relies solely on the AOA spectrum from a single station, making it highly practical for engineering applications.